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UNCLASSIFIED  ABSTRACT 

(U)  This  technical  report  describes  analyses  and  techniques  used  in  the  design 
and  evaluation  of  advanced  decoy  concepts.  The  work  described  addresses  both  the 
design  of  specific  penetration  aid  elements  and  the  formulation  of  techniques  for 
their  evaluation.  The  three  major  technical  areas  covered  in  this  report  are: 

1.  Investigation  of  a  penetration  aid  technique  that  degrades  the  measurement 
capability  of  the  radar  sense r. 

2.  The  design  of  a  computer  program  to  solve  the  decoy  design  problem  with 
flexibility  in  the  selection  of  optimization  criteria  and  constraints. 

3.  Studies  of  the  use  of  certain  discrimination  techniques  for  a  hard  point 

*  defense  system. 

This  appendix  to  this  report  contains  detailed  description  of  the  mathematical  and 
engineering  formulation  of  the  optimum  decoy  design  program. 
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l.o  I  INTRODUCTION  AND  PURPOSE  OF  PROGRAM 


1.1  PURPOSE  AND  GENERAL  REMARKS 

The  Optimum  Decoy  Design  Program  selects  optimum  decoy  configurations  which  meet 
specified  performance,  weight,  and  geometric  constraints.  Optimum  typically  means 
minimum  weight  for  a  given  set  of  performance  and  geometric  constraints  or  mini¬ 
mum  value  of  a  performance  parameter  for  a  given  weight  and  set  of  geometric 
constraints.  Mutually  exclusive  constraints,  for  which  there  is  no  acceptable 
solution,  are  identified  by  the  program.  The  program  is  capable  of  operating  in 
a  mode  which  allows  it  to  search  for  any  decoy  which  meets  all  the  constraints, 
(.solution-finder)  as  well  as  the  mode  where  tire  search  is  for  the  optimum  decoy 
within  all  the  constraints.  The  selection  processes  involve  the  minimization  of 
non-linear  functions  subject  to  non-linear  constraints.  The  analysis  portion  of 
the  program  includes  trajectory  calculations  (which  can  include  the  effects  of 
mass  loss,  nose  blunting,  thrust,  trailing  appendages,  and  angle  of  attack),  wake 
RCS  and  length  approximations,  and  effectiveness  operations.  The  program  is 
designed  primarily  for  sharp  and  blunt  cone  decoys  and  reference  reentry  vehicles. 
An  option  is  provided  which  allows  an  arbitrary  reference  reentry  vehicle  to  be 
utilized  provided  its  performance  histories  are  available  for  input. 

The  analysis  portion  of  this  program  does  not  represent  Avco's  most  accurate 
prediction  techniques  but  rather  represents  a  carefully  balanced  set  of  engineering 
approximations  which  are  utilized  to  maintain  the  required  compromise  between 
accuracy  and  running  time  for  configuration  selection  studies.  The  analysis 
portion  of  the  program  does  not  explicitly  contain  calculations  for  the  internal 
arrangements  of  the  equipment  inside  the  decoys  or  for  the  exoatmospheric 
observables.  The  effects  of  these  requirements  on  the  decoy  configuration  selec¬ 
tion  process  have  usually  been  included  as  geometric  constraints  on  the  decov. 

The  ADTECH  III  contract  resulted  in  a  computer  code  which  performed  the  solution¬ 
finding  operations  automatically  for  constraints  based  on  trajectory  matching 
requirements  (Ref.  III-l).  The  objective  of  the  ADTECH  IV  effort  in  this  area 
was  to  develop,  document,  and  deliver  a  computerized  decoy  configuration  selection 
program  which  can  select  optimum  decoys  considering  trajectory,  effectiveness, 
and  wake  performance  parameters.  This  objective  was  achieved  by  utilizing  and 
extending  the  results  obtained  in  the  ADTECH  III  computer  program.  The  ADTECH  III 
solution-finding  program  represents  the  first  steps  toward  the  application  of 
digital  machine  optimization  to  the  decoy  design  problem.  The  analysis  portion 
(function  generator)  of  that  program  contains  the  decoy  trajectory  calculations 
and  the  comparison  of  the  physical  characteristics  of  the  decoy  trajectory  with 
the  reference  reentry  vehicle  trajectory.  The  synthesis  portion  of  the  program 
utilizes  a  penalty-function  equation  to  transform  the  constrained  design  problem 
into  an  unconstrained  search  problem.  This  penalty-function  equation  was  formu¬ 
lated  to  be  compatible  with  the  later  extension  of  the  program  to  include 
optimization  capability.  Optimization  techniques  are  applied  to  the  transformed 
function  until  any  decoy  is  obtaineu  which  is  within  all  the  constraints.  During 
the  ADTECH  IV  effort,  the  ADTECH  III  program  was  extended  in  the  analysis  portion 
to  include  effectiveness  operations  and  observables  approximation  and  extended  in 
the  synthesis  portion  to  include  the  capability  for  finding  the  optimum  decoy 
design  which  is  also  within  all  the  constraints- 
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The  program  developed  during  ADTECH  IV,  when  given  an  initial  decoy  design,  the 
ranges  of  interest  of  the  decoy  parameters,  and  the  system  constraints  and  opt i 
mlzation  criteria,  will  find  the  best  vehicle  which  is  within  all  the  constraints 
or  will  establish  that  there  is  no  solution  to  the  problem.  For  example,  t  e 
program  is  able  to  find  the  minimum  weight  decoy  which  is  within  the  geometric 
and  observables  constraints  and  which  has  a  trajectory  that  matches  a  reference 
trajectory  within  specified  tolerances  (corridors).  The  program  also  is  able  to 
find  the  vehicle  which  minimizes  a  parameter  related  to  the  Avco  effectiveness 
model  (probability  of  discrimination)  and  is  also  within  specified  weight  an 
geometric  constraints. 

1.2  SUMMARY  OF  THE  COMPUTER  PROGRAM  FORMULATION 

The  major  elements  of  the  Optimum  Decoy  Design  Program  are  outlined  in  Figure  II1-1. 
The  kev  problem  in  formulating  this  design  program  has  been  the  integration  of 
many  subroutines  into  an  efficient  engineering  design  program.  These  subroutines 
fall  into  two  major  categories.  The  first  category  includes  those  subroutines 
which  perform  the  analysis  tasks  that  generate  vehicle  configurations  and  that 
generate  the  aerodynamic  and  trajectory  behavior  of  these  configurations.  In 
addition,  this  portion  determines  reentry  observables,  performs  a  comparison 
analysis  between  the  decoy  and  reentry  vehicle,  and  places  constraints  on  the 
observable  behavior  of  the  R/V  and  decoy.  These  constraints  are  generally  of  a 
corridor  type  or  of  an  effectiveness  parameter  type. 

The  general  decoy  conceptual  design  which  can  be  modeled  by  this  program  is  de¬ 
fined  by  the  thirteen  independent  design  variables  identified  in  Figure  HI-2. 
Besides  the  normal  blunt  cone  design  parameters,  provisions  are  made  to  cons id e 
variable  geometry  techniques,  in  which  the  decoy's  geometric  and  mass  charac- 
teristics^catTbe^changed^ln  bight  at  a  designated  altitude.  In  addlt  en  an  on 
generator  or  rocket  motor  can  be  incorporated  whcse  operational  characterise! 
are  controlled  by  the  analysis  portion  of  the  program.  The  decoy  design  problem 
then  is  the  search  for  the  numbers  defining  a  decoy  configuration  which  meets  all 
the  systems  requirements  and  constraints  and  is  optimum  with  respect  to  some 
specified  payoff  quantity.  This  analysis  portion  of  the  optimum  decoy  program  is 
thus  a  compendium  of  many  Avco  analysis  procedures  that  have  been  efficien  y 
combined  into  this  program. 

The  second  major  portion  of  the  optimum  decoy  design  program  is  the  synthesis 
section  which  employs  selected  optimization  techniques  to  find  optimum  s°lu^°ns 
to  the  decoy  design  problem.  The  optimization  techniques  chosen  for  inclusion 
°  this  program  L  Le  Davidon  Variable  Metric  Technique,  the  Rosenbrock  Rotating 
Coordinate  Technique,  and  the  one-  and  wo-v.rl.ble  Fibonacci. techn  q«s.  The 
formulation  of  the  program  will  be  discussed  in  more  detail  in  the  folio  g 

paragraphs . 

As  shown  in  Figure  III-l,  the  analysis  and  synthesis  portions  of  the  program  re¬ 
quire  inputs  which  define  the  problem  statement  and  the  options  required  or 
obtaining  a  solution. 

Ihe  inputs  to  the  program  Include  the  systems  requirements,  starting  conflgu- 

ration  and  the  ranges  of  allowable  parameter  variations.  A  solution  to  the 
problem  can  consist  of  any  decoy  configuration  which  meets  all  the  system  require¬ 
ments  and  constraints,  or  it  can  be  an  optlju.decoy  ten  gj  lo  whi  ^1^1 
all  the  constraints,  depending  on  whether  the  solution  find  g  p 

finding"  options  are  selected. 
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With  the  inputs  to  the  problem  specified,  the  trajectory  and  observables  of  the 
initial  decoy  are  calculated  in  the  analysis  portion  of  the  program  and  these 
data  are  compared  with  the  reference  reentry  vehicle  data  to  determine  the 
"performance"  of  this  particular  decoy.  The  effects  on  the  performance  as  a 
result  of  changing  the  parameters  of  the  decoy  configuration  slightly  are  also 
calculated  if  required  by  the  synthesis  portion  of  the  program,  lhese  effects 
ire  called  influence  coefficients  or  partial  derivatives  or  collectively,  the 
gradient  With  this  information  from  the  analysis  portion  of  the  program,  the 
synthesis  portion  of  the  program  changes  the  decoy  configuration  to  make  it  more 
acceptable  and  the  analysis  is  repeated  for  the  new  decoy.  This  process  is 
continued  until  a  decoy  is  determined  which  meets  all  the  stated  requirements  or 
until  it  is  determined  that  there  is  no  acceptable  decoy.  The  later  case  comes 
about  because  it  is  possible  to  state  a  problem  for  which  there  is  no  solution 
That  is,  the  problem  may  be  over-constrained.  The  program  is  designed  to  identify 
such  a  situation  and  indicate  that  no  solution  exists  for  the  stated  problem. 

In  addition,  the  program  defines  a  decoy  which  is  closest  to  being  acceptable. 

The  term  configuration  in  this  context  implies  the  specification  of  the  values  of 
the  design  variables.  An  example  Is  shown  in  Figure  II1-2  where  there  are  thirteen 
design  variables  which  must  be  determined.  This  decoy  initially  has  a  weight,  W,  , 
and  a  shape  defined  by  the  nose  radius,  RN]  ,  base  radius,  RBj  ,  and  length,  4. 

At  some  altitude,  ZTURN,  the  outer  shape  is  removed  and  the  decoy  immediately 
after  that  altitude  is  defined  by  the  parameters  W2 ,  R^,  R^,  and  L2.  In  addition, 

there  is  an  ion  generator  or  rocket  motor  which  initiates  at  some  altitude,  7vn, 
and  burns  out  at  some  altitude,  Zo(f  .  The  thrust  is  specified  at  some  value,  T. 
and  the  mass  loss  due  to  thrust  is  determined  from  the  thrust  and  the  specified 
value  of  specific  impulse,  I  .  The  decoy  design  problem  then  is  the  search  for 
the  numbers  defining  a  decoy  Configuration  which  meets  all  the  systems  require¬ 
ments  and  constraints. 

The  problems  Involved  with  the  specification  of  "suitable"  performance  or  systems 
requirements  is  not  within  the  scope  of  this  task.  The  systems  requirements  for 
input  to  the  program  must  be  provided  from  separate  studies.  For  example  velocity 
histories  for  a  reentry  vehicle  and  a  decoy  are  shown  in  Figure  1II-3.  The 
ouestion  of  whether  or  not  this  "match"  allows  the  decoy  to  perform  its  mission 
in  some  environment  is  considered  to  be  a  systems  problem  outside  the  scope  of 
this  computer  code.  However,  once  the  system  requirements  are  provided,  one  of 
the  essential  features  of  the  program  is  the  capability  for  determining  whether 
or  not  a  particular  decoy  meets  the  stated  requirements. 

In  addition  to  the  requirements  which  may  be  placed  on  the  trajectory  and  the 
observables  there  are  design  constraints  which  place  limits  on  the  weight  and 
geometry  of  the  decoy.  For  example,  the  length  and  base  radius  might  be  limite 
by  the  design  of  an  existing  deployment  system. 

The  search  for  an  acceptable  decoy  configuration  is  then  typically  bounded  by  a 
number  of  constraints,  some  trajectory  related,  some  observables  related,  and 
some  related  to  the  size,  shape,  and  mass  of  the  decoy. 

Provisions  are  included  to  constrain  each  of  the  design  parameters  between  lower 
and  upper  limits  if  desired.  The  program  can  accept  any  combination  of  the 
constraints  which  are  available  in  the  program. 
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SAMPLE  VELOCITY  HISTORIES 


The  above  sections  have  described  the  decoy  design  problem  and  the  general  pro¬ 
cedures  utilized  to  solve  the  problem.  In  the  following  paragraphs  the  actual 
formulation  of  the  computer  program  will  be  discussed  in  terms  of  its  modules  or 
subroutines,  starting  with  the  outermost  program  and  proceeding  in  to  the  actual 
trajectory  calculations.  The  main  program  (Figure  III-4)  accepts  the  input 
quantities,  allows  the  reentry  vehicle  data  to  be  input  or  calls  the  basic  function 
evaluation  subroutine  (F123)  the  reentry  vehicle  data,  allows  the  selection  of 
the  desired  search  subroutines  which  in  turn  call  the  function  evaluator  (FEV) 
and  provide  for  the  final  output. 

The  function  evaluator  (Figure  111-5)  screens  the  coordinates  (decoy  design  param¬ 
eters)  selected  by  the  optimizer  to  determine  if  they  are  compatible  with  the 
capabilities  of  the  trajectory  subroutines.  If  the  coordinates  are  not  acceptable, 
a  dummy  function  is  defined  to  encourage  the  optimizer  to  select  more  reasonable 
designs.  If  the  coordinates  are  acceptable,  the  matching  subroutine  (F123)  is 
used  to  compare  the  reentry  vehicle  data  with  decoy  data  calculated  using  the 
trajectory  and  observables  subroutines.  The  matching  subroutine  calculates  the 
corridor  fi  notions  illustrated  by  the  shaded  area  in  Figure  III-6.  The  corridor 
function  immediately  indicates  whether  or  not  the  decoy  velocity  history  is  within 
the  allowable  corridor.  If  the  function  is  zero,  the  velocity  aspects  of  the 
trajectory  are  acceptable;  if  the  function  is  positive,  the  velocity  history  is 
not  acceptable.  This  function  also  provides  a  measure  of  the  "amount"  the 
velocity  history  is  out  of  the  corridor.  Provisions  are  made  for  calculating 
corridor  functions  for  nine  observable  histories: 

a.  Velocity, 

b.  Deceleration, 

c.  Ballistic  Coefficient, 

d.  Wake  length  at  the  first  frequency, 

e.  Wake  length  at  the  second  frequency, 

f.  Wake  length  at  the  third  frequency, 

g.  Wake  RCS  at  the  first  frequency, 

h.  Wake  RCS  at  the  second  frequency, 

i.  Wake  RCS  at  the  third  frequency. 

In  addition,  nine  "effectiveness  integrals"  are  also  optionally  calculated  in  the 
matching  subroutine  for  use  in  the  probability  of  discrimination  calculation  in 
subroutine  EFFECT. 

The  formulation  of  the  miscellaneous  subroutine  (MISC)  includes  calculations  of 
the  average  density  of  the  decoy  which  may  be  constrained  to  aid  in  obtaining 
reasonable  packaging  situations.  This  average  density  calculation  is  a  first  step 
toward  the  possible  requirement  for  internal  packaging  logic  which  can  be  added 
into  this  subroutine  at  some  future  time.  Also  provisions  are  made  in  this  sub¬ 
routine  to  allow  the  parameters  describing  the  decoy  after  a  discontinuous  shape 
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change  to  be  constrained  relative  to  the  parameters  describing  the  decov  .just 
before  the  shape  change.  For  example,  the  vehicle  must  not  be  allowed  to  gain 
weight  at  the  shape  change. 

The  subroutine  called  EFFECT  contains  the  operations  for  defining'a  probability 
of  discrimination  based  on  any  desired  combination  of  the  nine  performance 
functions.  This  probability  of  discrimination  can  then  be  constrained  or  minimized, 
it  desired. 

The  calculation  of  the  penalty  function  equation  is  performed  in  the  function 
evaluator.  The  penalty  function  equation  provides  the  fundamental  link  between  the 
function  evaluator  calculations  and  the  synthesis  operations.  The  penalty  function 
equation  has  been  developed  to  handle  all  the  constraints  simultaneously.  This 
equation  is  an  adaption  of  the  work  of  Schmidt  and  Fox  (Ref.  I II — 2 ) .  An  Avco 
developed  programming  technique  for  this  equation  allows  any  combination  of 
quantities  which  are  available  in  the  program  to  be  constrainted . 

If  the  decoy  design  is  within  all  the  constraints,  the  penalty  function  will  be 
zero;  if  not,  the  penalty  function  will  have  some  positive  value.  This  formula¬ 
tion  immediately  provides  identification  of  acceptable  designs.  Information  about 
the  penalty  function,  !■'  ,  and  its  behavior  as  the  decoy  design  is  changed  is 
generated  and  used  to  proceed  from  a  starting  configuration  to  an  acceptable 
solution  and  then  to  an  optimum  solution,  if  desired. 

Typically,  the  behavior  of  the  penalty  function  is  sampled  for  each  of  the  design 
variables  being  considered  and  all  variables  are  then  changed  simultaneously  to 
proceed  to  a  new  configuration. 

An  illustration  of  the  use  of  the  penalty  function  is  shown  in  Figure  T.II-7  for 
two  independent  variables.  Conceptually,  an  optimization  technique  is  applied  to 
the  penalty  function  to  find  decoy  designs  with  smaller  values  of  the  function 
until  a  design  is  found  where  the  function  is  zero.  This  implies  that  the  final 
decoy  is  within  the  trajectory  and  observable  corridors  and  is  within  all  stated 
geometric  and  other  constraints.  This  process  is  repeated  with  a  sequence  of 
tighter  constraints  in  order  to  obtain  an  optimum  solution. 

If  a  minimum  is  located  where  the  function  is  not  zero,  this  implies  that  the 
stated  problem  is  overconstrained.  If  multimodel  functions  arc  suspected,  it  may 
be  advisable  to  start  the  search  from  several  different  starting  points  to  build 
confidence  that  there  is  not  a  better  solution  in  some  other  part  of  the  allowable 
search  volume. 

The  use  of  this  form  of  penalty  function  has  provided  an  immediate  identification 
of  acceptable  solutions,  a  procedure  for  identifying  impossible  problems,  and  has 
mathematically  transformed  a  constrained  optimization  problem  into  a  sequence  of 
unconstrained  optimization  problems. 

The  synthesis  portion  of  the  program  is  designed  to  perform  a  sequence  of  searches 
for  decoys  which  have  zero  values  of  the  penalty  function.  As  shown  in  Figure  III —4 , 
four  techniques  have  been  mechanized  for  selection  by  the  user  of  the  code  de¬ 
pending  on  the  nature  of  the  decoy  design  problem  and  on  the  experience  the  user 
lias  with  optimization  techniques.  The  search  techniques  consist  of  the  Davidon 
Variable  Metric  Technique,  the  Rosenbrock  Rotating  Coordinate  Technique,  and  the 
one-  and  two-variable  Fibonacci  techniques. 
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i’ll  is  optimum  decoy  design  program  has  been  utilized  on  several  programs  to  design 
decoys.  Its  utility  has  been  demonstrated  on  programs  like  the  Nike-X  Targets 
Program  and  it  is  being  utilized  in  examining  Tethered  Radar  Reflector  (TRR)  de¬ 
signs  in  the  program  of  that  name.  As  a  result  of  the  effort  pursued  under  this 
task,  the  penetration  aids  community  now  has  a  tool  that  economically  performs 
a  series  of  parametric  evaluations  that  normally  would  take  an  engineer  many 
months  of  effort.  Unlike  most  parametric  studies,  this  program  rigorously  chooses 
an  optimum  design  and  ensures  that  no  configuration  has  been  left  out  of  considera¬ 
tion  for  a  given  requirement. 

1.3  IDENTIFICATION  OF  DOCUMENTATION  ELEMENTS 

The  documentation  for  this  code  is  in  three  parts.  The  volume  entitled, 

"Numerical  Description  of  the  Optimum  Decoy  Design  Program"  (Appendix  I) 
contains  descriptions  of  each  subroutine,  descriptions  of  the  numerical  methods 
employed,  correlation  of  program  segments  with  their  function  descriptions  and  a 
complete  listing  of  the  source  program  and  preset  deck.  The  numerical  description 
volume  is  primarily  intended  for  the  professional  programmer  who  is  required 
to  understand  the  internal  detailed  operation  of  the  code.  The  second  part  is 
entitled  "User's  Manual  for  the  Optimum  Decoy  Design  Program"  (Appendix  II). 

This  manual  contains  detailed  descriptions  of  the  input  symbols,  input  sheets, 
sample  inputs,  and  corresponding  sample  outputs,  along  with  descriptions  of  the 
output  and  tests  of  the  correctness  of  the  operation  of  the  program.  The  user's 
manual  is  primarily  intended  for  immediate  reference  by  the  personnel  preparing 
inputs  for  the  code.  The  third  part  is  entitled  "Mathematical  Formulation  of 
the  Optimum  Decoy  Design  Program"  (Appendix  III)  and  contains  the  basic  equations 
and  procedures  which  are  mechanized  in  the  code.  This  part  is  of  use  to  all 
personnel  involved  with  the  code;  however,  it  is  primarily  intended  to  provide  a 
description  of  the  state-of-the-art  of  the  technology  utilized  in  the  code  which 
will  guide  the  current  users  and  will  serve  as  a  baseline  for  futher  improvement 
in  the  code.  The  three  parts  of  the  description  overlap  in  some  areas  in  order 
to  provide  more  complete  explanations;  however,  it  is  expected  that  a  careful 
user  of  the  code  would  read  all  three  parts  to  reduce  the  chance  of  misunder¬ 
standing  the  description  in  any  one  part.  In  particular,  the  User's  Manual 
presupposed  an  understanding  of  this  Mathematical  Formulation  document.  These 
three  parts  of  the  program  description  are  designed  to  be  independent  of  the 
lask-by-task  discussion  of  the  portion  of  the  effort  accomplished  during  ADTECH  IV 
which  is  included  in  the  main  body  of  the  ADTECH  IV  final  report. 

The  following  sections  of  this  volume  describe  the  formulation  of  the  Optimum 
Decoy  Design  Program  in  detail.  The  material  in  Sections  2.0  through  6.0  in  this 
volume  and  in  Sections  2.0  through  6.0  of  the  numerical  description  volume 
(Appendix  1)  are  organized  in  roughly  the  same  order.  The  breakdowns  of  the 
individual  sections  have  been  made  different  in  some  cases  to  best  serve  the 
needs  of  each  volume. 
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2.0  OPTIMIZATION  TECHNIQUE 


The  techniques  described  in  this  section  include  the  definition  and  utilization 
of  a  penalty  function  transformation,  the  calculation  of  the  gradient  of  the 
penalty  function,  and  the  description  of  the  various  search  techniques  which 
can  be  applied  to  the  penalty  function. 

2.1  PENALTY  FUNCTION  TRANSFORMATION 

The  penalty  function  used  in  this  program  allows  the  constrained  optimization 
problem  to  be  solved  as  a  sequence  of  unconstrained  problems.  This  section 
contains  descriptions  of  the  penalty  function  equation,  the  technique  for 
defining  a  dummy  function  if  the  current  search  point  is  outside  the  capability 
of  the  trajectory  modules,  and  the  techniques  for  tightening  a  constraint  after 
each  unconstrained  search. 

2.1.1  Function  Evaluation 

The  analysis  portion  of  the  program  (including  the  calculation  of  the  corridor 
functions,  miscellaneous  quantities,  and  effectiveness  parameters)  characterizes 
each  decoy  by  a  set  of  numbers  called  y;  .  The  selection  of  the  desired  com¬ 
bination  of  these  quantities  to  be  used  in  the  penalty  function  is  made  by  the 
user  by  means  of  input  codes.  The  penalty  function  equation  is  shown  in  Fig¬ 
ure  III-5 .  The  notation  used  will  be  defined  in  paragraph  2.1.3  of  this  volume. 
An  alternate  way  of  describing  the  calculation  of  the  value  of  the  penalty 
function,  F,  is  as  follows.  Let  F  be  the  sum  of  the  NP  individual  constraint 
terms,  f;  : 

NP 

i  =  1 

where  each  term  is: 

i  0.0  if  the  constrained  quantity  is  within  the  limits 


f  k,  ( |  amount  out  of  the  iimits  |  ) 1  if  not  within  limits 

where  k;  and  I  are  inputs.  If  the  penalty  function,  F,  is  zero,  the  decoy  is 
within  all  the  constraints.  If  the  penalty  function  is  not  zero,  then  its 
magnitude  represents  the  "amount"  that  this  decoy  is  unacceptable.  The  penalty 
function  thus  reduces  the  entire  set  of  performance  histories  and  design  con¬ 
straint  situations  to  a  single  number  (F).  Optimization  techniques  can  then  be 
applied  to  find  the  decoys  which  result  in  smaller  values  of  the  penalty  function 
until  a  decoy  is  found  which  has  a  penalty  function  value  of  zero.  This  would 
then  complete  one  "Solution-finding"  operation.  The  detail  of  the  process  of 
utilizing  a  sequence  of  solution-finding  operations  to  obtain  an  optimum  solu¬ 
tion  will  be  discussed  in  paragraph  2.1.3. 
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It  is  possible  for  an  optimizer  to  specify  that  a  decoy  be  evaluated  which  is  not 
physically  meaningful  or  which  is  outside  the  capability  of  the  trajectory  module  . 
h  Sp=e  has  been  set  bp  which  "screens"  out  such  situations  and  de- 

fines  a  value  for  the  penalty  function,  F,  which  is  equal  to  the  sum  of  the  last 
calculated  value  from  the  penalty  function  equation,!)  ,  to  the  amount  the  decoy 

failed  the  screening  test: 


Si|<LLi-Xi>2 


<X;  -  UL;>' 


where  S'  is  input-  and  the  parameters,  X;  ,  the  lower  limits  LL;  ,  and  the  upper 
limits'  UL-  are’ given  in  Table  III-I.  It  is  assumed  in  this  formulation  that 
the  reference  reentry  vehicle  and  the  initial  decoy  are  input  so  that  both  pass 
the  screening  test.  Otherwise,  it  is  likely  that  the  search  process  will  g 

unstable. 

2.1.3  Sequential  Optimization 

.  t- iniip  for  formulating  the  solution-finding  problem  was  discussed  in  para- 
*  £13“  S  This  technique  has  been  extended  to  allow  per  ornance 

and  design  constraints  to  be  included  in  the  search  for  an  optimum  decoy  config 
uration  &  The  extended  technique  effectively  transforms  the  constrained  optim 

ssssss 

following  form: 


F  =  kt  <Vl  -  UBp1 


E  i'i' 


+  <Vi  -  U<Y 


where: 

F  is  the  penalty  function  which  is  searched  for  a  zero  or  a  minimum  by 
the  optimizer 

kj  is  the  multiplier  (scale  factor)  for  the  quantity  being  minimized 
y  is  the  quantity  being  minimized 

UBj  is  an  upper  bound  on  the  quantity  being  minimized 

E  is  the  exponent  which  is  normally  2  but  can  be  set  to  1  for  special 
problems 

M  is  the  total  number  of  terms  in  the  penalty  equation 
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TABLE  HI-1 

SUMMARY  OF  SCREENING  LIMITS 


Conditions 

Lower  Limit 

Parameter 

Upper  Limit 

Units 

0.0 

RN, 

10000.0 

inches 

If  ZTURN  '-0.0 

0.0 

rN2 

10000.0 

inches 

0.0 

rNj  rb, 

0.6 

— 

If  ZTURN  >0.0 

0.0 

rn2/rb2 

0.6 

— 

3.0 

LAI 

168.0 

inches 

If  ZTURN  0.0 

3.0 

LA2 

168.0 

inches 

4.0 

T11FTA1 

40.0 

degrees 

If  ZTURN  >0.0 

4.0 

THF.TA2 

40.0 

degrees 

-10.0 

ZTURN 

ZO 

feet 

I  f  ZTU  RN  >0.0 

0.0 

W2 

W1 

pounds 

If  NTH  RUST  4  0 

0.0 

ISP 

10000.0 

seconds 

If  NTH  RUST =  I 

Z0FF 

ZON 

ZO 

feet 

If  NTHRUST  =  2 

TO 

TON 

TOFF 

seconds 

If  NTHRUST =  1 

ZST 

Z0FF 

ZON 

feet 

If  NTHRUST =  2 

T0N 

TOFF 

TST 

. 

seconds 
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ki  is  the  multiplier  (scale  factor)  for  each  constraint 

LHj  Is  the  lower  bound  for  each  constraint 

y |  is  the  value  of  each  item  being  constrained 

UH;  is  the  upper  bound  for  each  constraint 

vA>  is  |  A  if  A  is  positive 

I  0.0  if  A  is  zero  or  negative 

Note  that  the  user  should  input  a  lower  bound  for  the  quantity  being  minimized 
which  is  well  below  the  expected  minimum.  This  drops  the  first  "lower  bound" 
term  from  the  equation  and  thus  this  term  is  not  shown  in  the  equation  above. 

There  are  two  cases  to  be  considered.  The  first  is  where  the  quantity  being 
minimized  is  a  calculated  value  or  a  design  variable  and  the  second  is  where 
the  quantity  is  an  input  quantity  which  is  not  an  active  design  variable. 

For  the  cases  where  the  quantity  being  minimized,  yj ,  is  a  calculated  value  or  a 
design  variable,  the  procedure  is  to  input  the  exponents,  the  upper  and  lower 
bounds,  and  the  multipliers.  The  first  upper  bound,  Ubj  ,  is  selected  to  be 
reasonably  large  (based  on  experience  and  judgment).  The  multipliers,  kj  ,  are 
allowed  for  numerical  purposes  and  do  not  theoretically  enter  into  the  process 
at  all.  The  first  solution-finding  operation  is  initiated  from  an  input  initial 
decoy  configuration.  The  optimizer  modifies  the  configuration  until  every  con¬ 
straint  is  satisfied.  At  this  point  every  term  of  the  penalty  equation  is  zero 
and  the  value  of  the  penalty  function,  F,  is  therefore  zero.  The  parameters 
describing  this  acceptable  decoy  are  carried  over  to  describe  the  initial  con¬ 
figuration  for  the  next  solution-finding  step.  The  first  upper  bound,  UBj  ,  is 
reduced  for  the  second  solution-finding  operation.  Each  time  that  it  is  possible 
to  find  an  acceptable  decoy,  the  constraint  is  tightened  by  a  factor,  and 

the  process  is  repeated:  RI” 

<ubi>m  =  wrf  (y*)i 

After  a  series  of  solution-finding  problems,  the  constraint  will  have  become  so 
tight  that  a  solution  is  no  longer  possible.  In  this  situation,  the  search 
procedure  finds  that  F  has  a  minimum  and  that  it  is  not  possible  ti  find  a 
decoy  where  F  equals  zero,  thus  it  is  not  possible  to  find  a  decoy  which  satis¬ 
fies  all  the  imposed  constraints.  The  next-to-the-last  decoy  then  is  declared 
the  optimum  (within  the  tolerance  implied  by  the  factor,  WRF).  The  process 
described  here  is  the  one  mechanized  in  the  program;  however,  it  is  recognized 
that  it  may  be  desirable  in  the  future  to  incorporate  more  sophisticated  methods 
for  finding  the  smallest  compatible  value  of  the  upper  bound  UBj .  Approaching 
the  optimum  from  the  "acceptable"  side  has  been  found  to  be  a  practical  technique 
for  design  work. 

A  parallel  procedure  is  used  when  the  quantity  being  minimized,  y,  ,  is  an  input 
which  has  not  been  identified  as  a  design  variable.  After  each  successful 
solution-finding  operation,  the  value  of  the  input  Itself  is  changed: 
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\gain,  L he  solution-finding  process  is  repeated  until  it  is  no  longer  possible  to 
liml  n  solution.  The  next-to-the-last  decoy  is  then  declared  the  optimum  (within 
the  implied  tolerance). 

When  the  general  multivariable  optimizers  are  used,  the  exponents  in  the  penalty 
equation  must  have  a  value  of  two  so  that  the  first  derivatives  of  the  function 
with  respect  to  the  design  variables  will  be  continuous.  However,  since  the 
Fibonacci  techniques  do  not  have  any  continuity  restrictions,  the  exponents  can 
he  set  to  a  value  of  one  and  the  constrained  optimum  can  he  searched  for  directly. 
For  one-  and  two-var.able  problems  where  Fibonacci  techniques  arc  applicable, 
setting  the  exponents  to  one  and  providing  suitable  lj  multipliers  eliminates 
the  need  for  the  series  of  solution-finding  operations  discussed  above. 

2.2  GRADIENT  OF  THE  PENALTY  FUNCTION 

The  gradient  of  the  penalty  function  is  calculated  using  finite  differences.  The 
increments,  \X;  ,  are  inputs  to  the  program.  Each  element  of  the  gradient  is 
calculated  as: 


F(Xj  *  \X;)  -  F  (Xj) 
\X, 


2.3  SEARCH  LOGIC 

Four  separate  search  techniques  are  mechanized  for  possible  selection  by  the  user. 
The  first  is  the  Davidon  Variable  Metric  Method  for  Minimization  which  is  a 
multivariable  unconstrained  technique  which  utilizes  the  gradient  and  an  approxi¬ 
mation  of  the  matrix  of  second  partial  derivatives  to  locate  the  minimum.  The 
second  technique  is  the  Rosenbrock  Rotating  Coordinate  Method  which  locates  the 
unconstrained  minimum  of  a  multivariable  function  without  the  direct  calculation 
of  the  gradient.  The  third  technique  is  the  one-variable  Fibonacci  method  which 
locates  the  minimum  of  an  unimodel  function  within  a  specified  interval.  The 
fourth  technique  is  the  two-variable  Fibonacci  method  which  locates  the  minimum 
within  a  specified  region  by  obtaining  the  minimum  of  a  series  of  one-variable 
solutions.  These  four  search  techniques  will  be  discussed  in  more  detail  in  the 
following  sections. 

2.3. 1  Davidon  Method 

The  Variable  Metric  Method  for  Minimization  is  a  multivariable  technique  which 
uses  special  gradient  methods  to  locate  unconstrained  minimums  or  zeros  of  a 
function.  The  report  written  by  William  C.  Davidon  which  describes  this  method 
(Ref.  Ill —3)  is  reproduced  as  Appendix  1 1 1  —  1  to  this  volume  since  the  report  Is 
not  now  readily  available  (Ref.  TTI-4,  p.  20)  and  since  the  material  contained 
in  it  is  essential  to  the  understanding  of  the  method.  Additional  summary 
comments  and  notes  regarding  modifications  are  presented  in  the  following  para¬ 
graphs.  A  simplified  flow  chart  of  the  Davidon  Method  is  shown  in  Figure  III-8. 
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Figure  111-8  SIMPLIFIED  DAVIDON  FLOW  CHART 


This  method  is  a  modified  gradient  technique  having  more  sophistication  than  the 
first-order  steepest-descent  methods  and  far  less  computational  requirements  per 
slop  than  the  second-order  gradient  methods.  The  program  requires  an  initial 
starting  coordinate  (initial  design  point),  X  ,  an  estimate  of  the  inverse  of 
the  matrix  of  second  partial  derivatives  of  the  function  with  respect  to  the 
design  variables,  H  ,  and  a  procedure  for  evaluating  the  function, F  ,  and  its 
gradient,  G.  The  Davidon  method  selects  a  step  size,  A  ,  and  a  search  direction 
based  on  the  modified  gradient,  HG.  The  next  decoy  design  to  be  evaluated,  X  , 
is  determined  from  the  matrix  equation: 

X  X1’  -  A  HG 

The  step  size  is  estimated  by  the  program  (based  on  the  value  of  the  function, 
the  gradient,  and  the  estimate  of  H  to  overstep  or  bracket  the  minimum  along  the 
search  direction.  If  the  minimum  has  been  bracketed,  an  interpolation  is  made 
ter  the  approximate  minimum.  Based  on  the  behavior  of  the  function  during  the 
overstep  and  the  interpolation,  the  metric  H  is  modified  and  the  process  is  re¬ 
peated  from  the  best  point  available.  The  identity  matrix  is  typically  used  for 
the  initial  estimate  of  H  unless  other  information  is  available. 

it  is  the  use  of  first-order  calculations  (gradient)  to  improve  the  estimate  of 
the  second-order  parameter  H  which  gives  the  Davidon  technique  its  power  and 
efficiency.  If  the  penalty  function  becomes  zero,  an  acceptable  decoy  design 
has  been  obtained  and  the  problem  is  complete.  If  the  transformer-gradient 
becomes  less  than  a  tolerance,  -  ,  a  nonzero  minimum  has  been  located.  This  im¬ 
plies  that  it  is  not  possible  to  reduce  the  function  further. 

If  the  estimated  step  size  does  not  bracket  the  minimum,  the  Davidon  program 
modifies  the  H  matrix  and  selects  a  new  search  direction.  An  alternate  approach 
which  is  more  consistent  with  the  basic  theory  is  to  change  the  estimated  step 
size  until  the  new  point  does  bracket  the  minimum.  Although  this  approach  is 
theoretically  sound,  actual  experience  in  running  both  approaches  has  indicated 
that  the  original  approach  is  faster  for  most  problems. 

Another  variation  away  from  the  basic  theory  is  included  in  the  program  at  the 
point  where  the  minimum  in  the  search  direction  has  been  bracketed.  Calculations 
are  made  to  estimate  whether  the  function  is  behaving  in  such  a  way  that  a 
"perpendicular"  step  would  be  more  advantageous  than  back-tracking  to  inter¬ 
polate  for  the  minimum  in  the  search  direction.  This  possibility  is  attractive 
because  of  the  characteristics  of  certain  functions.  Safeguards  are  included  in 
the  program  so  that  the  process  is  stable.  For  example,  if  the  perpendicular 
step  is  worse  than  expected,  the  program  returns  to  the  normal  logic  and  inter¬ 
polates  for  the  minimum  in  the  search  direction. 

The  program,  as  mechanized,  differs  from  the  description  in  the  appendix  with 
regard  to  the  limits  on  the  maximum  size  of  the  step  length,  A.  In  the  program 
the  step  length  is: 

(2-° 

A  smaller  of  / 

I  -  M  <!/*«) 

where  M  is  an  input  which  may  be  adjusted  by  the  use»\ 
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me  random  step  operations  described  for  sebroutln.  STUFF  have  not  been. 
since  they  would  tend  to  interfere  with  the  sequential  constrained  opt  1ml  . 

process. 

2.3.2  Kosenbrock  Method 

The  Rosenbrock  Rotating  Coordinate  Minimization  Technique  (Reference  II 1-5)  is  an 

ilgorithm  Cor  selecting  trial  values  for  the  input  parameters  of  a  given  system 
model  in  such  a  way  that  a  function  of  the  performance  of  the  system  optimize  . 

fhe  technique  is  referred  to  as  "rotating  coordinates"  because  of  the  fashion  in 
which  the  variables  are  perturbed.  Rather  than  varying  the  Jnput  parameters  of 
the  system  model  one-at-a-time ,  this  method  rotates  the  coordinate  system  m 
parameter  space  so  that  one  axis  points  in  the  "best"  direction  of  search  The 
remaining  axes  for  exploration,  which  are  mutually  orthogonal,  are  obtained  with 
”  “schmldf  ortho  penalization  procedure.  A  search  Is  . :to» 
these  direction,  ooc  at  a  time  u,io6  the  lostc  which  s  toe  bed  below  to 
then  a  new  set  of  axes  are  developed.  A  geometric  interpretation 
systems  is  shown  in  Figure  II 1-9  and  flow  charts  for  the  technique  are  shown  in 

Figures  III— 10  and  III-ll. 

The  search  technique  involved  consists  of  taking  trial  steps  along  each  of  the 
coordinate  axes.  ’the  trial  Is  a  "success"  if  the  ^tieoal  -loe  * 
the  value  on  the  previous  trial;  otherwise  it  is  a  failur  .  p  . 

n-e  determined  in  the  following  manner.  An  initial  step  size  (AX)  is  an 
quantity^  After  a  successful  trial,  the  length  of  the  previous  step  1b  mu  Hi died 
by  a  constant  «,  («  >  D  and  this  is  added  to  the  previous  value  used  to  locate 
the  next  trial.  After  a  "failure,"  the  previous  step  size  is  multiplied  by  /  , 

(!)  1  <  0  and  this  is  added  to  the  previous  value  If  this  stepping  procedure 

produces  a  functional  value  within  a  specified  tolerance  level  (TOL .)  of  the 
lllTus  functional  value  the  step  is  called  a  "success  "  The  trials^in  a  g  .en 
direction  are  complete  when  there  has  been  a  "success  followed  by  a  failure. 

The  initial  step  in  the  new  search  (i.e.,  after  the  rotation  of  the  coordinate 
axes)  is  dependent  upon  the  total  of  the  successful  steps  from  the  previous 
rch  In  particular,  if  dn  is  the  algebraic  sum  of  all  successful  trials  in 

dlrectlon^then  to  W.1  trie!  of  the  next  search  in  the  e*  Irect  on 
wiLl  be  yd  (>  '•  0),  where  y  is  a  preset  constant.  Although  all  of  the  above 
constants,  « ,  ft,  V,  TOL,  are  preset,  they  may  be  input  as  different  values  to 
suit  the  need  of  a  particular  problem. 

After  a  set  of  trials  has  been  completed  in  one  direction,  the  program  searches 
lone  the  next  orthogonal  direction  until  all  N  directions  have  been  treated.  A 
new  set  of  directions  if  then  calculated.  All  of  the  trials  along  the  „ 

tions  and  the  subsequent,  calculation  of  a  new  set  of  directions  is  called  stage. 

The  rotating  coordinate  axes  are  related  to  the  parameter  axes  by  the  direction 

ostno  matrix  lQ„  I  ,  with  which  steps  in  a  given  direction  can  be  resolved  in 
cosine  matrix  (n  ,  matrix  so  that  each  step, 

parameter  changes.  For  the  first  stage,  iu J  v  For  each 

c  corresponds  to  a  change  in  only  one  of  the  system  parameters  Xp.  For 
subsequent  stage,  a  new  direction  cosine  matrix  is  computed  using  the  Gram- 
Schmidt  procedure  as  follows: 
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Lot  f,°  ,  £2° . .  £N°be  the  set  o£  ortho8onal  unit  vectors  defining  the  direc¬ 

tions1  in  the  original  stage Suppose  that  dj  is  the  algebraic  sum  of  all  success 
ful  steps  in  the  direction  etc.  Then  define  the  set  of  vectors: 

A1  dl  f"i°  '  d2  £2°  4  '  •  '  •  4  dN  £n° 

A2  d2  f2°  4  •  ■  '  •  +  dN  4j° 


An  dN  ^N° 

The  orthogonal  unit  vectors  f/  ,  f2l,  ...,^for  the  next  stage  are  now  obtained 
using  the  following  vector  equations: 


nl  A1 

*7  V  V 


(  An  •  f  j1  ) 


The  new  direction  cosine  matrix_is  obtained  by  taking  the  transpose  of  the  matrix 
comprised  of  the  components  of  V  along  the  parameter  axes.  With  the  new  coordi¬ 
nate  system  defined,  the  search  is  repeated  in  each  of  the  new  directions  in  turn 
The  result  of  applying  these  equations  several  times  is  to  ensure  that  q  lies 
along  the  direction  of  fastest  advance,  f2  along  the  best  direction  which  can  be 
found  normal  to  and  so  on. 


»2  A2  "  '  51 

?2  "  Vi  "2 1 


N  -  1 


i  =  1 

■?N  '  BN  !  i  bN  ' 


The  stopping  logic  is  based  on  the  value  of  the  function.  For  a  value  of  the 
function  that  is  undefined,  an  error  message  is  printed  out  and  the  search  is 
stopped.  The  same  thing  happens  if  the  total  number  of  function  evaluations  for 
a  given  set  of  constraints  equals  an  input  limit.  If  the  value  of  the  function 
is  zero  or  if  three  succeeding  functional  values  are  within  a  defined  interval  of 
each  other,  a  solution  is  considered  to  have  been  found.  Considering  the  diagram 
below,  the  intervals  for  successive  functional  values  which  define  a  solution  are 


<Fs 


Fs  +  2)  <  DEL*  Fs  and 


F s  +  1  “  ^st2 

FS-FS+I 


<  RATU 


where  DEL  and  RATU  are  input  limits. 
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The  functional  values,  Fs  ,  Fs  +  j  ,  Fs +  2,  are  those  obtained  when  the  trial  steps 
have  been  completed  for  all  the  system  orientations. 

2.3.3  One-Variable  Fibonacci  Method 

The  Fibonacci  search  technique  is  a  search  scheme  for  finding  the  maximum  or 
minimum  of  a  one-variable  function  within  defined  limits.  The  function  has  to 
be  at  least  piecewise  continuous,  single  valued,  and  also  have  only  one  optimum 
(i.e.,  maximum  or  minimum)  within  the  interval.  These  restrictions  define  a 
unimodal  function  (pp.  10-13  of  Reference  III-6) .  The  initial  interval,  Lj  f  is 
defined  by  an  upper  bound,  B  ,  and  a  lower  bound,  A  ,  where 

Lj  =  B  -  A 

Either  the  number  of  function  evaluations,  F  ,  to  be  made  during  the  search  or 
an  end-of-search  accuracy  limit,  Ac  ,  defined  in  terms  of  a  number  of  independent 
variable  units  away  from  the  actual  maximum  or  minimum  within  the  interval  also 
has  to  be  given. 

The  technique  is  based  on  direct  comparison  of  values  of  the  function,  which  are 
used  to  exclude  parts  of  the  search  interval  (Figure  III-12) .  The  nlacement  and 
comparison  of  the  points  is  continued  until  the  interval  is  sufficiently  small 
or  until  the  function  is  zero.  The  points  ere  located  so  as  to  maximize  the 
interval  to  be  excluded  at  each  step.  The  three  cases  shown  in  Figure  111-12 
illustrate  the  comparison  and  exclusion  steps.  The  left  end  of  the  first  case 
can  be  excluded  since  the  minimum  must  be  either  in  the  center  or  right  sections. 
The  right  end  of  the  second  case  is  excluded.  The  third  case  shows  that  both 
ends  can  be  excluded  if  both  values  of  the  function  are  exactly  the  same.  Flow 
charts  of  this  technique  are  presented  in  Figures  111-13  and  III-14. 


u 
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The  location  of  the  first  two  points  in  the  initial  search  interval  is  dependent 
upon  t’.ie  total  number  of  function  evaluations  to  be  made  during  the  search.  If 
the  accuracy  is  given,  a  trial  Fibonacci  number,  F^-p  ,  is  defined 


The  actual  Fibonacci  number  is  obtained  through  an  iteration  process  where 


F  ■>  2 


|;n  f'n  -  1  *  l;n  -  2  "  3'  ‘4-  '  * 

The  iteration  process  is  continued  until  the  first  value  of  n  at  which  Fp  Fn-j- 
is  found.  The  number  of  function  evaluations  to  be  used  within  a  given  accuracy 
limit  is  then: 

Nf  n  -  2 

Likewise,  if  the  number  of  function  evaluations  is  given  the  associated  Fibonacci 
number  is  found  by  applying  the  above  process  in  reverse  order. 

The  two  evaluation  points  in  the  first  interval  are  located  as  follows: 


and  the  corresponding  functional  values,  Yj  and  V2  are  computed.  Comparing  the 
two  functional  values  and  choosing  one  of  the  associated  evaluation  points  to  be 
a  new  end  point  yields  a  new  interval  for  computing  the  next  evaluation  point. 

For  instance,  in  the  search  for  a  maximum  the  evaluation  point  associated  with 
the  smallest  functional  value  becomes  an  end  point  for  the  new  interval  containing 
the  other  evaluation  point  and  the  remaining  end  point.  Then  a  new  evaluation 
point  is  placed  in  the  new  interval  a  distance  of  Ax  units  from  the  end  point 
retained  from  the  previous  interval. 


i  =  3,  A,  -  *  -  ,  Nj^ 


where  the  subscript  j  is  the  number  of  the  particular  evaluation  point  being 
found  for  that  interval.  The  new  functional  value  for  the  interval  is  computed 
and  compared  with  the  functional  value  retained  from  the  previous  interval  so 
that  a  new  interval  can  be  determined.  This  process  is  repeated  until  the 
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appropriate  number  of  function  evaluations  have  been  made.  maximum  or  m 

mum  ot  the  function  is  obtained  by  directly  comparing  the  functional  vaiue  of 
the  final  evaluation  point,  the  evaluation  point  retained  from  the  last  *  , 

and  the  new  end  point  for  the  last  interval. 

2.3.4  Two-Variable  Fibonacci  Method 

Tilt’  Fibonacci  scorch  technique  for  o  two-varlablc  function  In  a  bounded 
i,  the  application  of  the  one-varlablc  Fibonacci  search  o  each  variable ,  1 

one  of  the  variables  to  be  the  secondary  Independent  variable i  and  the  other  to 
be  the  primary  independent  variable,  the  method  employed  in  this  technique  is 
the^opt  I  mlz.it  ion  of’the  primary  independent  v.rlab.e  durln*  enc h  .  ep  o 

tsrss-jLS 

primary  independent  variable  is  optimized.  Then,  the  secondary  value  and 

£  33  zzzr 

variable.  The  table  is  then  searched  for  the  optimum  secondary  value,  producing 
also  the  associated  optimum  primary  value. 

The  restriction  of  unimodallty  (see  description  of  one  variable  Fibonacci)  applies 
III  a  stricter  sense  for  the  two  variable  Fibonacci  technique  Not  only  does  the 
two  variable  function  have  to  be  unlmodal,  but  the  associated  one-variable  func- 

independent  variable  is  important.  For  the  case  of  a  unimodal  two  variable  tu 
tion  where  both  choices  for  the  secondary  independent  variable  lead  to  multi 
modal  Sn'variabL  functions,  the  two-variable  Fibonacci  search  technique  may 

not  be  reliable. 
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3.0  BASIC  ANALYSIS  CALCULATIONS 


The  analysis  calculations  involve  the  basic  calculations  required  for  a  single 
vehicle  and  the  comparison  and  effectiveness  calculations  which  are  required  to 
evaluate  the  degree  of  simulation  of  two  vehicles.  The  basic  calculations  are 
discussed  in  Section  3.0  and  include  the  trajectory,  wake,  and  miscellaneous  cal¬ 
culations.  The  comparisons  are  discussed  in  Section  4.0  and  the  effectiveness 
model  operations  are  discussed  in  Section  5.0.  The  analysis  calculations  are 
performed  in  order  to  define  the  quantities,  y ; ,  in  the  penalty  function  equation 
of  Section  2.1. 

3.1  TRAJECTORY  CALCULATIONS 


The  trajectory  calculations*  are  incorporated  into  the  ADTECH  Decoy  Design  program 
to  provide  the  flight  characteristics  of  the  candidate  decoy  designs.  Three 
methods  of  computing  the  trajectory  are  available  (a)  particle  trajectory,  (b) 
rotational  three  degree  of  freedom  method,  and  (c)  simplified  a  gle  of  attack 
solution.  The  effects  of  ablation,  thrust,  and  angle  of  attack  oscillation  can 
be  included  in  the  trajectory  calculations.  The  program  operates  within  the 
following  geometric  and  flow  constraints: 


Cone  half-angle 
Body  length 
Surface  temperature 
Bluntness  ratio 
Altitude 

Frce-stream  Mach  number 
Angle  of  attack  at  300  kft 
Angle  of  attack  below  150  kft 


4°  <  0  <  40° 

C  — 

0.25'  <  La  <  14’ 

1000°  R  £  T  <  6000°  R 
0.0  £  A  £  0.6 
0.0'  £  Z  <  400,000' 

5.0  <  <30.0 

0°  <  a  <_  20° 

0°  <  a  <_  0„ 


The  mathematical  description  of  the  trajectory  program  is  presented  in  this  part 
of  the  report.  The  description  is  broken  up  into  sections,  each  containing  an 
individual  area  of  analysis  involved  in  the  program.  Since  there  is  no  direct 
correspondence  between  the  sections  in  the  mathematical  description  and  the  sub¬ 
routines  in  the  numerical  description,  the  brief  description  of  each  section  given 
below  lists  the  subroutines  that  contain  the  area  of  analysis  being  described. 


The  preliminary  calculations  section  (3.1.1)  contains  the  flow  field,  geometry, 
and  thrusting  calculations.  This  includes  the  free-stream,  edge,  and  wall  flow 
field  properties,  geometry  before  and  after  shape  change,  and  all  of  the  thrusting 
parameters.  The  subroutine  invclved  is  PRELIM. 


The  heating  section  (3.1.2)  contains  the  calculation  of  the  aerodynamic  heating 
rates  at  various  stations  along  the  body.  The  subroutine  involved  is  AERODY . 

•Developed  for  the  U.  S.  Army  Material  Command,  Redstone  Arsenal,  Alabama,  under  Contract  l)A-01-021-AMC-9003d(Y). 
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The  ablation  effects  (3.1.3)  section  contains  the  calculation  of  the  sidewall  and 
nose  recession  rates,  mass  loss  rates  at  various  stations  along  the  body,  and  the 
total  vehicle  mass  Joss  rate.  The  subroutines  involved  are  MASSLO,  EVIL,  NOSEBL, 
and  T0MAL0. 

The  angle  of  attack  section  (3.1. A)  contains  the  calculation  of  the  oscillatory 
angle  of  attack  of  the  vehicle.  The  subroutine  involved  is  a  portion  of  ROTATE. 

The  drag  section  (3.1.5)  contains  tiie  calculation  of  the  flight  drag  coefficient. 
The  effects  of  ablation  and  angle  of  attack  are  included  in  the  calculations. 

The  subroutine  involved  is  DRAGCO . 

The  equations  of  motion  section  (3.1.6)  contains  the  calculation  of  the  derivatives 
of  the  flight  trajectory  parameters  and  the  integration  of  the  equations  of  motion 
to  determine  the  flight  trajectory.  The  subroutines  involved  are  DEREQ,  TEQUAT, 
ADMARK,  and  a  portion  of  ROTATE. 

The  values  of  the  curve-fit  coefficients  called  \  and  B;  in  the  equations  of 
Section  3.0  are  defined  in  the  source  listing  of  subroutine  ZPRS  in  Appendix  1. 


3,1.1  Preliminary  Calculations 

The  preliminary  calculation  of  the  geometric,  flow  field,  and  thrusting  parameters 
is  done  in  the  subroutine  PRELIM.  The  operations  involved  can  be  broken  up  into 
seven  majo:  areas  of  calculation:  (a)  body  geometry  definition,  (b)  aerodynamic 
coefficients,  (c)  atmospheric  free-stream  properties,  (d)  wind  tunnel  free-stream 
properties,  (e)  flow  field  properties,  (f)  thrusting  quantities,  and  (g)  weight 
increments.  In  the  following  analysis  each  of  the  seven  areas  will  be  presented 
separately,  along  with  a  listing  of  the  inputs  necessary  to  perform  the  calculations. 

To  define  the  geometric  properties  of  the  body,  the  following  inputs  are  necessary: 
,cone  half-angle,  0C  (degrees),  nose  radius,  Rfj(feet),  and  base  radius,  Rg  (feet). 
The  sharp  cone  slant  length  is  given  by: 


sin  H 


(ft) 


the  axial  length  from  the  stagnation  point  by: 

r  rn  1 

LA  -  L  co»flc  -  —  (l-*in0c)J  (ft) 

and  the  surface  distance  from  the  nose  of  a  blunt  cone  by: 


'B 


sin 


e . 


kn_ 

tan  0. 


(ft) 


"C  -  c 

where  0C  has  the  units  of  radians.  The  bluntness  ratio  is  defined  as. 
rN 


rb 

The  base  diameter  is  D  «  2Rfi  (ft) 
and  the  reference  area  is: 


'ref 


n  R, 


(ft.2) 
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The  inputs  necessary  to  compute  the  aerodynamic  coefficients  are  the  cone  half- 
angle,  0C  (degrees),  nose  radius,  Rjg  (feet),  base  radius,  Rr  (feet) ,  bluntness 
ratio,  '1  ,  and  the  nondimenslonal  center  of  gravity  location  from  the  nose,  Xc,g.  D* 
The  center  of  gravity  location  is  input  in  tabular  form,  as  a  function  of  altitude, 
with  the  appropriate  value  being  used  for  the  altitude  point  under  considerate n . 
The  partial  derivative  of  the  normal  force  coefficient  with  respect  to  angle  of 
attack  is  given  by: 


Cm  =  1.92  cos  20c  i 
« 


—  (A  cos  ftr)2 
2  c 


•p|, „  center  of  gravity  location  is  determined  by  multiplying  the  base  diameter  by 
the  appropriate  value  of  the  nondimenslonal  center  of  gravity  location  from  the 
input  table: 


The  partial  derivative  of  the  pitching  moment  coefficient  with  respect  to  angle 
of  attack  is  given  by: 


^  2[  1  -  (A  cos  tfc)?]  -  3  A  cos3  Gc  i  1  -  ( A  cos  0C)2] 

3  tan 

1  i 

—  (  A  cos  0  ) “ 

2  _ 

,  r  ♦ 

-  A  sin  0(  cos2  [  1  -  (A  cos  0C ) 1 1  ^ 

Using  the  above  results,  the  nondimenslonal  center  of  pressure  distance  from  the 
nose  of  the  body  becomes: 


(RN"Xc.R.)eos2  °c 


'CP 


M«  Xc.g. 


where  the  appropriate  value  for  - —  is  taken  from  the  input  table 


Tha  partial  derivative  of  the  pitching  moment  coefficient  with  respect  to  the 
pitching  rate  can  either  be  input  or  calculated  as  follows: 
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2 


A  co  s  ’  () 


1  -  ( A  cos  0C ) 2  ]  - 


4  A  cos'  (I 


3  sin2  t'l 


I  1  -  ( A  cos  0  ) '  ] 


1  -  (A  cos  0)4 


2  sin1"  0_ 


A  cos  9, 


(^)J 


2  C  X  -  Rn(  1  -  sin  0C)  I  , 

- - 2[  I  -  ( A  cos  0  ]  -  3  A  cos*  0  l  1  -  (A  cos  <L)2  I  r 

ARBtan(9c  I  c  c  ’ 

/  co*  rtc  \  2  | 

♦  --  )  l  1  -  <  A  cos  )  2  1  l  Vc  R  -  RN  (  l  -  sin  0C  ) )  2  ^ 

To  conduct  a  flight  analysis  on  a  body,  the  atmospheric  free-stream  properties 
must  be  computed.  The  necessary  inputs  are  the  altitude,  Z  (feet)  and  the  free- 
stream  velocity,  (ft/sec).  For  a  nonstandard  atmosphere,  a  table  of  free- 
stream  density,  (slug/ft3)  and  speed  of  sound,  a^  (ft/sec)  as  functions  of 
altitude  must  be  input  to  provide  these  two  values  at  the  altitude  point  being 
considered.  For  a  standard  atmosphere,  the  appropriate  values  of  the  free-stream 
density  and  the  speed  of  sound  are  obtained  from  the  1962  standard  atmosphere 
(subroutine  ARFDT2) .  Changing  units  on  the  £  ee-stream  density  to  lbm/ft3  gives: 

P*,  - 

Free-Stream  Temperature  (°  R) : 

Tk  »  4.16  x  1(T4  a  J 

2 

Free-Stream  Viscosity  (lb-sec/ft4): 

(32.2)  (2.27  *  10“8)(Tx)1-5 
_  _____ 

Free-Stream  Reynolds  Number  based  on  the  Axial  Length  of  the  Body: 

P-.j  u~  LA 


Free-Stream  Pressure  (lb/ft2): 


P*  P-j  RT» 
Free-Stream  Mach  Number: 
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To  perform  the  drag  calculations  for  wind  tunnel  conditions,  the  free-stream  Mach 
number,  M^,  free-stream  Reynolds  number  per  inch,  Re^Cl/in),  and  the  free-stream 
total  pressure,  P(  (lb/ft2),  are  necessary  inputs.  Changing  the  units  on  the 
free-stream  Reynolds  number  to  (1/ft),  one  has: 

Re^f  -  12.0  Rc^ 

The  Free-Stream  Pressure  (lb/ft2): 


1  +  M  2 


y/y-  1 


The  Free-Stream  Temperature  (  R) : 


1 


2*°  Re~f  |  1.490825  x  10-8  R 
4(  198.6)  Re„ 


P  N1 

QO  OC 


f  V  1.490825  x  10-8  R 
The  Free-Stream  Density  (lbm/ft^): 


1.490825  x  10  R 
1/2  | 

\ 


'~1  RT^ 

3 

or  changing  Units  on  the  Free-Stream  Density  to  slug/ft  : 

Pacj 

PK  ~  T2.174 

The  Free-Stream  Speed  of  Sound  (ft/sec): 

v'tT 


0.020)96 


The  Free-Stream  Viscosity  (lb-sec/f t^) 


2'i . 


= 


(32.2)  (2.27  x  10“8)  (T^  ) J*5 
198.6 


and  the  Free-Stream  Velocity: 


*=  aw 


The  flow  field  properties  at  the  edge  of  the  boundary  layer  are  determined  with 
the  aid  of  the  free-stream  and  geometric  properties  calculated  above.  Ihe  cal¬ 
culations  are  as  follows: 


The  Stagnation  Enthalpy  (Btu/lbm) 

Hs  =  0.24  t  (2.0  x  lO-'MlC 

2 

Stagnation  Pressure  (lb/ft  ): 


PS  =  p~ 


(y  *  i) 


-  1  f(l  -  y)  +  2  y 


2 

Free-Stream  Dynamic  Pressure  (lb/ft  ): 


=  ~  P' c 


Frce-Stream  Reynolds  Number  based  on  the  Sharp-Cone  Slant  Length: 


P«j  L 


Edge  Velocity  (ft/sec): 


Ue  [l  -  (M^sinfl,)1-9 


Edge  Pressure  (lb/ft  ): 


i  2.5  +  8-M^  sin  6C  ( 

p  =  P  1  +  2.8  (M  sin  0  —  — —  '  T“ 

e  *  i  00  c  [l.e  16.M,*,  sin  f/c  J 

Edge  Temperature  (°  R) : 


1  4-  0.0966  M„  sin  0  4-  0.2267  ( sin  0C)‘ 


for  M„  sin  9  <  5.7 

OO  C 


T  =  T  X 

C  OC 


4  2 

EE 

i=o  k=o 


A (90  +  i  4-  5k  sin  0C1 


•  r  p~  i 

i  _ 

.2116.. 


for  sin  &c  >  5.7 
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Edge  Density  (lbm/ft^): 


H  T 


Edge  Mach  Number : 
U_ 


2x. 


Edge  Viscosity  (lb-sec/f t") 

(32.2)  (2.27  x  10-8)  (T.) 


1.5 


c  Te  +  198.6 

and  Edge  Reynolds  Number  based  on  the  Surface  Distance  from  the  Nose  of  a 
Blunt  Cone: 


Pe  Ue  S 


Re. 


If  the  transition  altitude,  ZTR  ,  is  not  input,  it  is  calculated  in  the  following 
manner: 


1 


1 


£  £  £  A (300  +  j  +  3i  +  6k)  (0C)'1  (L)i  (A)k 


i  =  o  i  =  o  k=  o 


for  0  <  15°  and  A  <  0.3 


JTR 


1  2 


A(300  +j +  3i  +  6k)(15.)‘(L)i  (0.3)k 


i=o  i=o  k=o 


for  6  >  15°  and  A  >  0.3 


\ 

In  the  following  calculations  the  appropriate  value  of  wall  temperature  is  used: 


Tw(2,  8)  for  Z  ^  ZTR 
Tw(3 ,  8)  for  Z  <  ZyR 
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For  the  above  wall  temperatures,  the  notation  is  as  follows: 

Tw(2,  8)  =  wall  temperature  at  the  maximum  diameter  point  of  a  sharp  cone  in 

laminar  flow. 

Tw(3,8)  =  waii  temperature  at  the  maximum  diameter  point  of  a  sharp  cone  in 

turbulent  flow. 

These  two  values  are  calculated  in  the  ablation  effects  portion  of  this  write  up. 


Constant  Pressure  Specific  Heat  at  the  Edge  of  the  Boundary  Layer  (Btu/lbm0  R) 


0.2398  for  Tc  <  "00°  R 


te 


A(  105  +  i)  <Te)‘  for  700°  R  <  Te  <  5000°  R 


A(lll)  +  A  (1 12)  Te  for  Te  >  5000°  R 


Constant  Pressure  Specific  Heat  at  the  Wall  (Btu/lbm  R) : 


0.2398  for  Tw(j,  8)  <  700°  R 
5 


Y  A ( 105  +  i ) [ Tw ( j  ,  8)]  1  for  700°  R  <  Tw(j,8)  <  ‘5000°  R 
i  =  o 

A  ( 1 1 1)  +  A  (112)  Tw(j,  8)  for  Tw(j,8)  >  5000°  R 
Dimensionless  Wall  Enthalpy: 


C 


Pw  Tw(j,  8) 
33-86 


Wall  Viscosity  (lbf-sec/f t^) 


2*. 


Pw 


(32.2)  (2.27  x  10-8 )  [ Tw ( j ,  8)1  1,5 
T_7j,8)+  198.6 


Chapman-Rubensin  Constant  for  Edge  Conditions: 

Pw  Te 
Pe  Tw(j,8) 
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Viscous  Interaction  Parameter: 


\ 


Free-Stream  Chapman-Rubensin  Constant: 


/'  w  V 

C, 

c  Tw(j,  8) 


Hypersonic  Rarefaction  Parameter: 


\1 


[TJ7 

V  rhs 


The  calculation  of  the  thrusting  parameters  is  initiated  by  the  input  of  a  thrust 
history  table.  This  table  may  take  one  of  two  forms:  (a)  thrust  as  a  function  of 
altitude,  or  (b)  thrust  as  a  function  of  time.  For  the  convenience  of  running  a 
number  of  cases  with  the  same  profile  but  at  different  thrust  levels,  the  table 
is  in  terms  of  a  nondimensional  thrust,  Th/THq.  The  dimensional  thrust  along  the 

trajectory  is  then  obtained  by  multiplying  the  nondimensional  thrust  table  values 
by  an  input  reference  thrust  level,  Tjj  .  Also  required  as  input  are  the  thrust 

offset  angles  and  distances.  These  are  the  angular  misalignment  of  the  thrust 
vector  in  the  yaw  direction  ( <]>a  )  and  in  the  pitch  direction  (0  ),  measured  in 
that  order  in  the  body— fixed  coordinate  system,  and  the  distance  the  thrust  vector 
is  offset  from  the  two  coordinate  axes  that  are  normal  to  the  body  centerline 
(  \Z  and  \Y).  The  thrust  at  specified  altitudes  or  times  in  the  trajectory  is 
given  by: 


TABLF 


(LBF) 


where  linear  interpolation  is  used  to  determine  the  proper  value  for  altitudes  or 
times  that  are  between  the  tabulated  nondimensional  thrust  values.  The  effective 
thrust  acting  on  the  vehicle  is  then: 


TH  "  1HCW  -  Ae  (LBF) 

The  components  of  the  effective  thrust  in  the  body  coordinate  system  are: 

Tji  =  Tjjr  COS  0(1  COS  \ji 

THV  TH  cos  °a  sin  ‘A 

TH7  TH  sin  °a 
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The  moments  generated  by  thrust  offset  are  given  by: 


MX  =  -  AY  Th?  -  AZ  THy 
my  'z  thx  •  <xc.M.  -  La)THz 


MZ  *  <Xc.g.  -  LA>THY  -  AY  THX 

The  instantaneous  total  weight  of  the  vehicle  is  determined  by  subtracting  the 
weight  loss  increments  due  to  the  altitude  or  time  step  along  the  trajectory  from 
the  total  weight  of  the  vehicle  determined  in  the  previous  step.  For  the  case  of 
the  first  step  along  the  trajectory,  the  total  weight  used  is  the  initial  weight 
of  the  vehicle.  The  necessary  inputs  are  the  total  vehicle  weight  from  the  pre¬ 
vious  step  (  W0  ) ,  the  weight  loss  due  to  thrusting  for  the  present  step  (AWTH), 
and  the  weight  loss  due  to  ablation  for  the  present  step  (AW).  The  instantaneous 
vehicle  weight  at  the  end  of  any  trajectory  step  is  given  by: 

W  =  WQ  -  AWth  -  AW  (lb  ) 

and  the  associated  vehicle  mass  is: 

W 


3.1.2  Heating 

The  subroutine  AERODY  calculates  the  heating  rates  at  different  stations  along  the 
body  surface  in  continuum  flow.  (See  para.  3.1.5.)  For  a  sharp  cone  the  body 
stations  are  the  stagnation  point  and  the  maximum  diameter  point.  For  a  blunt 
cone  the  body  stations  are  the  stagnation  point,  tangent  point,  20  percent  station, 
40  percent  station,  60  percent  station,  75  percent  station,  90  percent  station, 
the  maximum  diameter  point  for  laminar  and  turbulent  flow,  and  the  sonic  point  for 
turbulent  flow  only.  The  percentage  stations  are  located  according  to  the  initial 
axial  length  of  the  cone.  The  notation  used  to  identify  the  different  heating 


values  is: 

Qd,  1) 

=  Stagnation  point  heating 

Q(2,  i) 

=  Laminar  heating 

Q (3 ,  i) 

=  Turbulent  heating 

Q(4,  1) 

=  Sonic  point  heating  in  turbulent  flow 

where 

i  =  i 

Implies  tangent  point  on  a  blunt  cone 
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j  =  2-7  Implies  the  20  percent,  40  percent,  60  percent,  75  percent, 

90  percent  stations  and  maximum  diameter  point  on  a  blunt  cone 

j  =  8  Implies  the  maximum  diameter  point  of  a  sharp  cone 

All  of  the  heating  rates  calculated  in  the  section  have  the  units  Btu/ft2-sec. 
The  quantities  needed  to  calculate  the  heating  rates  are  the  f roe-stream  density, 
(lbm/ft3),  f  ree-stream  Mach  number,  M^,  free-stream  velocity,  (J<v,  (ft/sec), 
cone  angle,  0C  (degrees),  bluntness  ratio,  A,  nose  radius,  ]<iy(feet),  cone  axial 
length,  (feet),  free-stream/stagnation  pressure  ratio,  P^/Ps ,  dimensionless 
stagnation  enthalpy,  HS/RT0  ,  stagnation  pressure,  Ps  (lb/ft2),  sharp-cone  slant 
length,  L  (feet),  and  transition  altitude,  7tr  (feet).  All  the  above  quantities 
are  calculated  in  the  subroutine  PRELIM. 


The  stagnation  point  heating  for  sharp  and  blunt  cones  is: 


1.76  x  104 


31 


U 


3.15 


.  0.002375 


,  2.6  1(T 


1/2 


for  A  =  0.0 


Qd,  1) 


1.76  x  104 


3  15 


10.002375  RN  /  \2.6  x  104 


1/2 


for  A  >  0.0 


For  a  sharp  cone,  the  edge-stagnation  pressure  ratio  at  the  maximum  diameter  is 
given  by: 


0.0331  exp  [  0.0064  0C-0. 33  (M^~  5.0)  °-85  ] 


+  4.68  x 


10-4  (0c)  1-88032 


If  Z  >  Z-pg,  the  laminar  flow  heating  rate  for  a  sharp  cone  at  the  maximum  diameter 
point  is  calculated  as: 


Q(2,  8) 


r  0.9736 

0.5142  [Logc(Hs/RT0)] 

[  10] 

0.9664  +  0C  |  5.28  x  10~3  +  2.88  x  10“4  0C] 
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For  Z  <  ZTR,  the  turbulent  flow  heating  rate  for  a  sharp  cone  at  the  maximum 
diameter  point  is  calculated  as: 

Kj  =  0.9  +  0.02  0,. 

K,  0.69  *  L A  1 3.18  x  10~2  -  6.9  x  10"4  LA1 

I'j'  “ 


Q(3 1  8)  - 


0.745  [  Loge  (HS/RT0)1 


0.6(K1t)(K3t) 


,0.8122 


For  calculating  the  appropriate  heating  rates  on  a  blunt  cone,  the 

sis  /sir  s  zrzfxrzz  -r  ° 

between  the  spherical  nose  and  the  conical  body  -  is  given  by. 

/v\  RN(l-sin0c) 


If  the  body  is  to  have  a  discontinuous  shape  change  at  some  point  in  the  trajectory 
(LwrJEed  iy  the  shape  change  altitude,  ZTmll).  the  following  Is  true  for 

7  s  ZyuRN  : 

La 


For  7  <  ZXURN  : 

la2 

\La  = 

LA 

where  L A !  is  the  axial  length  before  shape  change  and  is  the  axial  length 
after  shape  change.  The  body  stations  are  located  as  follows: 

20  percent  station 

=  1  -  0.8  ALa 


40  percent  station 


1  -  0.6  ALa 


60  percent  station 


I  -  0.4  ALa 
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75  percent  station 


1  -  0.25  ALa 


90  percent  station 
X 


L  i 


1-0.1  ala 


Maximum  diameter  station 


The  pressure  distribution  over  a  blunt  cone  is  determined  as  follows 


/  j  La  tan2  0c 

then  the  sharp  cone  maximum  diameter  pressure  is  used: 


Or  if 


1.13  Rn 

l.A  tan2  0C 

then  the  pressure  distribution  is: 


P 


e 


for  0C  <  20° 


E 

If  =  n 


A  (135  +  n  +  3i  +  9k) 

1.13  ~/k 

La  tan  20c  ^ 
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p. 


p« 

p. 


1  £.  i. 

n  -  0  i  0  k  =  0 


A  (41  +  n  f  2  i  +  6k) 


(0.174  ocy 


10.  \  1  \ 

—  1  0.2  Loge 


1.13  R 


N 


La  tan2  Bc 


-H 


for  0C  >  20° 


If  Z  >  ZTD,  the  laminar  heating  rate  is  calculated  for  the  blunt  cone.  At  the 


tangent  point,  the  blunt  cone  laminar  heating  rate  is: 

\  /p 

Q ( 2  ,  1)  =  1.732  Q ( 1  ,  1)  ^  0.007789  +  1.849 


1.6832 


+  0.841 


J  \ 


The  heating  rates  at  the  various  stations  along  the  conical  frustum  of  a  blunt 
cone  in  laminar  flow  are: 


0.5142  f  Loge(  HS/RT0)1 


0.9736 


Q(2,  j) 


[101 


[  0.9664  +  5.28  x  10-3  0C  +  2.88  x  10  4  02) 


LN  \ 

1  +  1.782  I  -  -  2.008 
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For  Z  ZTR  the  turbulent  heating  rate  is  calculated  for  the  blunt  cone.  The 
turbulent  flow  heating  rate  at  the  sonic  point  on  a  blunt  cone  is: 


Q(4,  1) 


3760.  (p^)0-8  dU3-45 
(0.002375)0'8  (10) 13,8  (Rn)0-2 

for  Z  <  115,000  ft. 

3760 (/)  )0.8[u  |  (2.254  +  2.246  x  10-5  Z  -  1.469  *  10” 10  Z2  t-  3.671  *  10-15  Z  3  ) 

(0.002375) 0,8  (Rn)°i2  [10|4(2-254  f  2.246  x  10-5  Z-  1.469  x  10~10  Z2  f  3.671  x  10-15  Z3) 

for  Z  >  115,000  ft. 


The  tangent  point  heating  rate  for  a  blunt  cone  in  turbulent  flow  is  given  by: 


0.8 


13460 


rnu- 


57.3 


°-2  \  2.375  xl0~3/ 


4.1 

6.0 


n. 


j_\  °-4 

6  / 


[  0.0001  UJ3'« 


for  Z  <  115,000  ft. 


Q (3,  1) 


13460 


P* 0, 


0.8 


rnu- 


°c 

57.3 


0,2  \2.375  x  10~3 


4.1 


'1- 


s/l 


!  N  0.4 

T 


[  0.0001  U, 


j  (2.254  +  2.246  x  10-5  Z  -  1.469  x  10-10  Z2  +  3.671  x  10-15  Z3) 


for  Z  >  115,000  ft. 
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The  heating  rates  af  the  various  stations  along  the  conical  frustum  of  a  blunt 
cone  in  turbulent  flow  are: 


K1T  =  0.9  +  0.02  0Q 


R 


N 


2T 


i  1.0  for  -  >  0.2 

j  RB 


1 0.6  ,  2.0 


for 


<  0.2 


K3T  =  0.69  +  0.03 18L  -  0.00069L2 


K5T  =  0.69  +  0.0318  La  -  0.00069  LA2 
K6t  =  0*6 
K.  7  t  ~  1.0 


KgT  =  0.69  +  0.0318 


rn  + 


-  0.00069 


LA 


RN  + 


sin  6C 


KgT  =  0.901  -  0.867 


1.13  Rn 
3  La  tan2  9C 


+  0.966 


1.13  Rn 

3La  tan2  0C 


III-48 


QO,  j>  =  -< 


1 21 16 


0.8  n  I  i  ,0.8122 

0.745  I  Log  (H  /RT  )] 

[10 


(K1t)  (K2t)  (K3t)  (K4t) 


(1 

V'-A 

\  113  % 
j  3LAtan2flc 

ii' 

0-8  /  ,0.8122 
0.745 1  Log  (H  /RT)] 

10 

V*/ 

/  2116 

1 

(kll)  (k2t)  (ksT)  (k9T) 

,o,  /r 

V  M3Rn 

V-A, 

i  LA  ta"2  °c 

i 

M 

°'8  fjo  0.745tLoge(Hs/RTo)]°-8122 

(K1t)  (K6t)  (K7t)  K8T) 

for 


1.13  R 


N 


3.1.3  Ablation  Effects 


lA  ta°2  °c 


”  ?  s  eels 

rate  is  then  obtained  by  integrating  the  individual  mass  loss  rates  over  the 
surface  of  the  vehicle.  The  steady  state  ablation  method  entails  an  iterative 

cSe°  s  rra°  eSim?i:arS 5tfqUati0nS  f0r  the  SUrfaCe  t-peraSS^hi  -S1^ 

ene^gy  considerlMon,  maSS  transfer  “te  S°lutl°n  includes  the  following 

losses  I  2  I  convective  energy,  conduction  flux,  surface  radiation 

sses,  and  sublimation  energy.  Nose  blunting  is  also  included  in  the  category 


III-49 


of  ablation  effects.  The  equations  involved  in  the  mass  loss  rate  calculations 
provide  reasonable  answers  only  for  normal  reentry  conditions.  The  purpose  of 
these  mass  transfer  rate  equations  is  to  obtain  a  reasonable  approximation  to  the 
effects  of  ablation  on  drag.  Consequently,  the  results  of  the  ablation  effects 
equations  should  not  be  used  in  heat  shield  design. 

The  subroutine  EVIL,  in  conjunction  with  the  subroutine  MASSLO,  is  used  to  cal¬ 
culate  the  mass  loss  rate,  m (1,  j)  (lbrn/f t2-sec) ,  surface  recession  rate  s(i,  j) 
(ft/sec),  and  the  wall  temperature  Tw(i,  j)  (°  R)  at  the  appropriate  body  stations. 
The  definitions  of  the  different  subscript  combinations  for  the  above  quantities 
are  given  in  the  drag  calculations  section  (3.1.5).  As  was  mentioned  above,  the 
method  employed  to  determine  the  mass  transfer  rate  is  heat  shield-material  de¬ 
pendent.  For  OTWR  in  turbulent  flow,  and  I.Ta  ,  Teflon,  or  any  other  input  material 
in  laminar  or  turbulent  flow  the  steady  state  ablation  method  is  used.  For  OTWR 
in  laminar  flow,  and  Carbon  Phenolic,  or  Phenolic  Nylon  in  laminar  or  turbulent 
flow,  curve  fits  for  m(i,  j)  and  s(i,  j)  as  functions  of  the  cold  wall  heating  and 
Tw(i,j)  as  a  function  of  s(i,  j)  are  utilized. 

To  perform  the  mass  loss  rate,  wall  recession  rate,  and  wall  temperature  calcula¬ 
tions  the  following  material  property  constants  are  needed: 

ft  J,  ft2 .  fly  ft  4  >  Href'  «•  P2'  F’  -  NSL  -  NGL’  NST>  NC,T>  CP2’  CPG’  and 


For  the  heat  shield  materials  mentioned  abov»  .  the  material  property  constants  are 
given  in  subroutine  CHNTBL.  Also  needed  is  the  stagnation  pressure,  Ps  (lb/ft  ), 
an  initial  guess  at  the  wall  temperature,  TWq  (°  R),  and  the  appropriate  heating 
rates,  Q(i,  j)  (Btu/ft2  sec). 

The  iterative  steady  state  ablation  method  proceeds  as  follows.  For  the  first 
pass  through  the  equations,  the  wall  temperature  is  set  equal  to  the  initial  guess 
at  the  wall  temperature:  Tw(i,  j)  =  TWq  .  Then,  for  each  subsequent  iteration,  the 
wall  temperature  calculated  at  the  end  of  the  previous  iteration  is  used 

25000. 

\T  =  0.005  Tw(i,  j)  +  ,r  ,, 

Tw0,  |) 


S(i,  j)  =  Pi  Tw(i,  j)  + 


Lo*el^2lTw('*  i)^3  |-04/Tw<i>  i> 


Loge(l.  11057  x  107)  -  1. 

1112  x  105/ Tw(i,  j) 

J 201.8340(K1)2  +  756.2732  (^-}j  Kj  -  12.7657  Kj 

2 

h  •  «  Gr,)] 
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Vj  13.654  f  469.585  K2 

y2  0.256012  •  0.005558  K2 

>',  5.345  x  10~6  -  4.27  x  10-'  K2 

The  nondimensional  wall  enthalpy  is  given  by: 

_  Y\  4  12  Tw(i’  v  2 

"w  33.86 


E,  - 

1  35.89 


-0.037  for  laminar  and  turbulent  flow  at  stagnation  point 

0,0  for  laminar  flow  at  all  other  body  station 

0.0  for  turbulent  flow  at  all  other  body  stations  or 
turbulent  sonic  point  flow 

0.0  for  laminar  and  turbulent  flow  at  stagnation  point 

-0.185  for  laminar  flow  at  all  other  body  stations 

-0.502  for  turbulent  flow  at  all  other  body  stations  or 
turbulent  sonic  point  flow 


Qof  Q(i,  i)(Ej)ei  (e2)'2 


m-51 


The  new  wall  temperature  to  be  used  in  the  next  iteration  is  then: 

AT  for  L  -  R  <  0  and  AT  >  0 

or  L  -  R  >  0  and  AT  <  0 

_0.5  At  for  L  -  R  <  o  and  At  <  0 

or  L  -  R  >  0  and  AT  >  0 

The  iterative  calculations  end  and  a  solution  has  been  found  once 

I  1.0 

|  L  -  R  |  <  <  or 

|  0.01  L 

or  if  100  iterations  have  been  made.  Then,  appropriate  wall  temperature  is 
the  Tw(i,  j)  used  in  the  last  iteration  and: 

rh(i,  j)  =  s(i,  j)(p2  +  M 
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Tin  curve  fit  solution  for  phenolic  nylon  heat  shield  material  in  laminar  or 
turbulent  flow  is  as  follows: 


10 


m(i,  j)  = 


j  l  -  2.5228  i  7.3759  *  10-3  Qfi,  j)  ] 
for  Q(i,  j)  £  100. 

-  1.62367642  x  10~3  +  1.78922793  x  10-4  Q(i,  j) 
i  1.32113696  x  10-8  [  Q(i,  j)  1  2  -  5.087475  x  10~12  [ Q(i,  j)  1 1 
for  100.  <  Q(i,  j)  £  3000. 


Q(i,  i) 


1.0  - 


1700 

ir 


1845.  +  11 


■  (i)J 


^  for  Q(i,  j)  >  3000. 

The  wall  recession  rate  for  phenolic  nylon  is  ratioed  from  the  wall  recession  rate 
of  OTWR  in  laminar  flow. 


s(i,  j) 


[  s ( i ,  i)(3TWR  ^  ^  i3 
m(i,  i )  OTWR 


Tw(i,  i) 


6944.74035  +  645.367146  (  Log,0  s(i,  j)  ] 


-  148.589173  [  LogjQ  s(i,  j) 


-4 


for  s(i,  j)  >  10 


9030.  +  1756  [  Log10  s(i,  j)  ] 


for  £(i,  j)  <  10 


-4 
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The  curve 


fit  solution  for  OTWR  heat  shield  material  in  laminar  flow  only  is 


0.0  for  Q(i,  j)  <  13.0 


-0.01929  +  0.00015  Q(i,  i) 


for  13.0  <  Q(i,  j)  <  15.6 


1.27424339  x  10-3  +  1.3607167  x  10”^  Q(i,  j) 


•1.09091516 


x  10-6  t  Q(i,  j)  ] 2  +  7.98275747  x  10  9  [  Q(i,  j)  1  3 


m(i>  j)  = 


-1.65210579  x  10-11  [Q(*.  ))  1  4 


for  15.6  <  Q (i,  j)  <  250. 


-  1.05650025  x  10-J  +  7.61118699  x  10  5  Q(>.  )) 


+  3.342517  x  10-8  [Q(i,  j)  1  2  -  6.91682422  x  10  12  [ Q(i,  i)  1  3 


for  250  <  Q(i,  j)  <  3000. 


3640  +  8.1 


for  <2(i,  j)  >  3000. 
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The  wall  recession  rate  Is: 


s(i,  i) 


-5.23741  x  10-4  +  1.6115  *  10-6  Q(i,  j) 
for  Q(i,  j)  <  1000. 

-  1.1119676  x  10-4  f  4.03376719  *  10-7  Q(i,  j) 

+  9.70131261  <  10"10(  0(i,  j)  l2  -  2.45527504  x  10“ 1 3  t  Q (i ,  j)  1 3 
for  1000  <.  0(i,  j)  <  30  00. 

m(i>  i) 

P2  +  A  p 

for  Q (i ,  j)  2  3000. 


and  the  wall  temperature  is  given  by: 

Tw(i,  j)  =  6346.34912  +  550.628796  [ Log  10  S(i,  j)  1 

+  19.6585366  (Log10S(i,  j)  ]  2 


For  carbon  phenolic  heat  shield  material  in  laminar  flow,  the  ablation  character¬ 
istics  are  given  by: 


m(i,  j)  - 


0.0  for  Q(i,  j)  <  8.0 

6.5788  x  10'4  -  1.59773  x  10-4  Q(i,  j)  +  9.8485  x  Iff6  [0(i,  j)  ] 2 

for  8.0  <  Q(i,  j)  <  11.3 
-  9.9516436  x  10-4  +  9.78022  x  10-5  Q(i,  j) 

for  11.3  <  Q(i,  j)  <  23.1 
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For  carbon  phenolic  heat  shield  material  in  turbulent  flow,  the  curve  fit  solution 
for  the  ablation  characteristics  incorporates  a  portion  of  the  first  iteration  from 
the  steady  state  ablation  method.  As  in  the  iterative  solution,  an  initial  trial 
value  is  chosen  for  the  wall  temperature:  Tw(i. j)  *  Tw  then,  the  solution  pro¬ 
ceeds  as  follows:  ° 


\T  =  0.005  T„  + 
wo 


25000. 


SpD  =  /?!  Tw^  + 


^/T*„ 


K,  - 


Loge(l.  11057  x  107 )  -  1.11  12  x  IoVt, 


201.8340(Kj)2  +  756.2732 


2117 


Kj  -  12.7657  Kj 


Kj  +  4 


2117 


>'j  =  13.654  +  469.585  K2 
)'2  =  0.256012  +  0.005558  K2 
y 3  =  5.345  x  10-6  -  4.27  x  10-7  K2 
The  nondimensional  wall  enthalpy  is  given  by: 


>'l  +  72  Twn  +  F3(Tw  )■ 


Hw  = 


33.86 


Ej  = 


35.89 


0.349 


.  RT„ 


T 


E,  = 


0.349 


,  RT„ 


+  17.945 


=  0.95  - 


H. 


’ref 


RT 
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k  -0.037  for  stagnation  point  in  turbulent  flow 

el  =  'i 

i  0.0  for  sonic  point  and  all  other  body  stations  in  turbulent 
'  flow 

^  0.0  for  turbulent  flow  at  the  stagnation  point 

e  B 

j  -0.502  for  sonic  point  and  all  other  body  stations  in  turbulent 
\  flow 

Qof  =  i)(Fi)Cl  (E2)'2 
SpD  Hs[p2NST  +  ApNGTl 

Oof 

«5b  =  e~<7  +  0-618  f2) 

Qt  =  4.7853  X  10“13  <(TW  )4 

Qc  =  E3  Qof 

qnet  =  0c  -  Qt 

Q*  =  Cp^  (6160  -  Tw  )  +  F  +  Nst 

Qnet 

ri,(i,  j)  =  __ 

.  ..  0.8  m(i,  j) 

s(i,  |)  =  - 

P2 

Tw(i,  j)  =  6687.41134  +  446.431845  (Log  1Q  S  (i,  j) 3 
+  12.0991623  [Log10S(i„  j)  3 2 
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The  subroutine  NOSEBL  is  used  to  calculate  the  derivatives  with  respect  to  time 
of  the  nose  radius,  base  radius,  and  side  walls  in  order  to  determine  the  shape 
change  for  a  body  in  continuum  flow  with  a  constant  cone  half-angle.  For  cases 
where  the  bluntness  ratio  is  1  ;ss  than  or  equal  to  10"3,  the  body  is  considered 
to  be  sharp.  In  the  following  calculations  the  appropriate  wall  recession  rate, 
based  upon  the  heat  shield  material  being  considered,  must  be  employed.  The  only 
input  necessary  for  these  calculations  is  the  wall  recession  rate,  S(i,  j)  which 
was  computed  previously.  For  laminar  or  turbulent  flow  and  sharp  or  blunt  cones, 
the  rate  of  change  of  the  nose  radius  is  given  by: 

S(l,  1)  sin  0Q 

Rkt  —  - - - - 

1  -  sin  0C 

For  a  sharp  cone,  the  maximum  diameter  sidewall  recession  rate  is: 

\  0.7071  S(2,  8)  for  Z  >  ZTR 
Sw  , 

J  0.8706  S(3,  8)  for  Z  <  ZJR 

For  a  blunt  cone,  the  maximum  diameter  sidewall  recession  rate  is: 


\  S(2,  7) 
j  S(3,  7) 


for  Z  >  Z-, 


for  Z  <  Z-, 


Then,  for  both  flow  conditions  and  both  body  configurations,  the  base  radius 
recession  rate  is  given  by: 


The  subroutine  TOMALO  is  used  to  determine  the  rate  of  change  in  weight  due  to 
ablation.  This  is  accomplished  by  integrating  the  mass  loss  rates  for  the 
different  body  stations  over  the  surface  of  the  body.  The  necessary  inputs  are 
the  mass  loss  rates,  m(i,  j)  at  the  various  body  stations,  the  axial  locations  of 


the  body  stations, 


£). 


,  and  the  axial  length  of  the  body,  LA  .  For  a  sharp 


cone,  the  total  mass  loss  rate  is: 


-  2.9618  v  2  (L.)2  -  m(2,  8) 

cos 

c 


for  Z  >  Z-, 


tan  6C 

3.04(2.0) 0,2  - —  (La ) 2  m (2,  8) 


for  Z  <  Z-, 
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For  a  blunt  cone,  the  total  weight  loss  rate  Is  composed  of  the  mass  loss  rate 
from  the  spherical  nose  cap  and  the  mass  loss  rate  from  the  conical  frustum.  The 
total  mass  loss  rate  on  the  spherical  nose  is: 


The  conical  frustum  contribution  to  the  total  mass  loss  rate  is: 
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3.1.4  Angle  of  Attack 


The  angle  of  attack,  which  is  incorporated  into  the  drag  calculations,  is  deter¬ 
mined  by  the  subroutine  ROTATE.  Two  methods  are  available  for  calculating  the 
angle  of  attack  a.  uncoupled  rotational  three-degree-of-freedom  method  and 
b.  simplified  angle  of  attack  method.  In  both  methods,  the  objective  is  to 
define  an  angle  of  attack  value  (a)  from  the  oscillatory  variations  in  angle  of 
attack  that  occur  from  perturbing  forces  and  moments.  For  the  three-degree-of- 
freedom  method  the  necessary  inputs  ire  the  pitch  and  yaw  Euler  angles,  0a  and 

u,  the  rolling  component  of  rotation,  P,  the  maximum,  a ',  and  minimum,  a"  , 
perturbed  angle  of  attack  values  from  the  previous  oscillation,  the  minimum,  tL  , 
and  maximum,  ,  time  limits  on  the  period  of  oscillation,  the  time  of  the  first 
maximum  in  angle  of  attack,  ^  ,  the  time  of  the  second  maximum  in  angle  of  attack, 
7 2  >  period  of  oscillation,  t'  ,  the  viscous  interaction  parameter,  y,  cone 

half-angle,  0C  ,  and  the  bluntness  ratio.  The  total  angle  of  attack  is  given  by: 

a  '  =  cos-*  [  cos  0Q  cos  i/f  ] 

The  angle  of  attack,  a  ,  forwarded  to  the  drag  calculations  is  dependent  upon  the 
flow  regime,  cycle  time,  rolling  rotation  rate.  For  the  strong  interaction 
regime  (x  >  *up) : 

a  =  a' 

For  laminar  and  turbulent  continuum  flow  (  y  <  yu^) : 

(a  '  for  t '  =  0 

a'  for  t '  >  tH 


—  for  t '  <  t,  and  P  =  0 

n  L 

=5 

/  **0 1 

a  +  a 

— - —  for  t '  <  tL  and  P  4  0 


For  the  case  where  the  period  of  oscillation,  t',  is  between  the  two  limiting 
cycle  times,  an  integrated  angle  of  attack  is  calculated.  As  noted  in  the  drag 
section,  this  angle  of  attack  is  used  only  to  modify  the  forebody  pressure  drag. 
Thus,  for  tL  <  f  <  tH  ,  the  angle  of  attack  is  calculated  as  follows: 

A2  =  [ Aj 2  +  Aj50c  +  Ajg  9C2]  +  [  A21  +  A24  9c  +  A27  9  2  ]A 
+  [A30  +  A 33  9C  +  A36  0c2]A2 


A3  =  t  A13  +  A16  ^c  +  a19  1  +  t  a22  +  A25  ^c  +  A28  1  A 

+  [A3J  +  A34<?c  +  Kyj  0c2jA2 
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Chen, 


(  A2  a  +  A3  a  ) 


The  simplified  angle  of  attack  method  is  a  straightforward  approach  requiring  no 
integration.  ThUngle  o£  attack  supplied  to  the  drag  calculations  is  dependent 
upon  the  period  of  oscillation,  t  ,  of  the  vehicle.  For  this  method  the  period 
is  calculated  below.  The  necessary  inputs  are  the  derivatives  of  the  pitching 
lent  efficient  with  ,  and  «  ,  and  c^.  derivative  of  nomal  force  coeffi- 

cient  with  a,  CN  ,  the  minimum  time  limit  on  the  period  of  oscillation,  tL  , 

initial  flight  path  angle,  yf0  ,  initial  Euler  angle  in  the  pitch  direction,  0a  , 
drao  coefficient,  Cn  ,  altitude  variation  of  the  axial  moment  of  inertia,  I  , 
vehiclc^nstantaneous  weight,  »,  and  the  geonettic  and  flow  field  properties. 

The  calculations  are  as  follows: 


/9Z  «  -Log^ 


(  e‘J-) 

\  0.076474/ 


T(  .  .  ...  1  altitude  (Zm  Z  )  is  being  considered  then  the  drag  coefficient 

£  1  0.8.  Otherwise  the  drag  coefficient  obtained 

from  the  drag  calculation  section  is  used. 


32.174  rr  Rg2  Z 
^K1  °  4  /3Z  W  1  sin  yf  f 


4RB2cMqW  ’ 
"  CNa  +  32.174  I  . 


2.174  *RB2Z2  r(CNn-CD)^|sinyfo 


VK2  “  ”  ; 

2W  (/3Z  sin  yf  f 


2RBffcMfl 

32.174  I 


The  period  of  oscillation  is  given  by: 


-2»rZ 


Pz  u»  sin  YL  VaK2 
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The  total  angle  of  attack  equations  necessitate  the  calculation  of  Bessel  func¬ 
tions  of  the  first  and  second  kind  of  the  zeroth  and  first  orders.  These  Bessel 
functions  are  to  be  functions  of  S  where: 

s’  -  2  vp~(aK1  +  aK2 ^ 

The  Bessel  functions  calculated  are: 

Y0  (S),  Y^S),  J0(S),  J;  (?) 


For  the  initial  altitude  (Z  =  Z0)  only,  the  following  calculations  are  made: 


'C2 


\ 


2.3769  x  10-3  ^  ^  AK1  J0  (S) 
PZ  +  Z0  eXP  (AK1  Poc') 


0a  pz  Ji  (s')  e*p(-  AK1  P*o) 


2  Z„ 


I 


\  Sfi Z  M  M  .  I 

( Yj  (S)  J0  (S)  -  J!  (S)  Y0(S)1  J 


-1 


2  Z„ 


*C1  = 


^  ^AK1-Ac2Yo(S) 
o 


Vs) 


The  angle  of  attack  ratio  at  the  initial  altitude  is: 
P»AK1 


shnfK 


K2  P« 


For  all  other  altitudes  the  angle  of  attack  ratio  is: 
Poo  AK1 


J7 


'  V  AK2  Poo 

The  total  angle  of  attack  is  then  given  by: 


[  Aci  Jo<S)  +  AC2  Yo(~)l  ePooAK1 


for  t  >  t. 


a'  =  (  0.63661  dn 


for  t  <  tj  and  Z  =  Z0 


0.6366  1  0 


(a/a0) 
“o  (a/aj ) 


for  t  <  tj^  and  Z  <  Z0 
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and 


n  a ' 


3.1.5  Drag 

The  analytical  expressions  for  the  drag  coefficients  are  ordered  according  to 
the  flight  regimes  encountered  during  the  reentry  process:  free  molecular 
transition  flow,  strong  interaction  flow,  laminar  continuum  flow,  and  turbulent 
contiruum  flow  respectively.  Since  the  theoretical  developments  in  each  regime 
are  based  on  different  analytical  models,  discontinuities  exist  at  the  regime 
interfaces.  To  provide  a  continuous  drag  history,  fairing  techniques  are  used 
between  the  flight  regimes.  The  viscous  interaction  parameter,  y,  the  rarefaction 
parameter,  yj,  and  the  transition  altitude ,  Zq-R  ,  dictate  which  flow  or  fairing 
regime  is  to  be  used  in  the  drag  analysis.  The  range  of  applicability  of  the 
given  body  parameters  and  flight  conditions  for  this  drag  model  are: 


•  Cone  Ha If- Angle 


4°  <  0C  <  21° 


Body  Length 
Surface  Temperature 
Bluntness  Ratio 
Altitude 

Free-Stream  Mach  Number 


1 '  <  La  <  14 ' 

1000°  R  <  Tw  <  6000°  R 

0.0  <  A  <  0.6 

0.0  '<  Z  <  400,000' 

5.0  <  <  30.0 


Angle  of  Attack  at  300  kft0°  <  a  <  20° 
Angle  of  Attack  <  150  kft  0°  <  a  <  $c 


Several  body  and  flow  parameters  are  needed  to  initiate  the  drag  calculations. 
They  are  cone  half-angle,  0C  (degrees),  angle  of  attack,  a  (degrees),  free-stream 
velocity,  U„  (ft/sec),  altitude  Z(ft.),  and  any  two  of  the  following:  bluntness 
ratio  A ,  base  radius  Rg(ft.),  or  nose  radius  Rfj(ft.).  Also  needed  are  the  wall 
temperature,  Tw  (i,j)  (°R),  mass  loss  rate.m  (i,j)  (lbm/ft^sec) ,  the  wall  en¬ 
thalpy,  Hw,  and  the  total  mass  loss  rate,  W(lbm/sec).  The  wall  temperature, 
wall  enthalpy,  and  mass  loss  rate  are  calculated  in  the  subroutine  EVIL. 


The  subscripts  i  and  j  refer  to  specific  flow  regimes  and  body  stations.  The 
subscript  combinations  used  in  the  drag  calculations  are  as  follows: 

1.1  -  stagnation  point 

2.1  -  tangent  point  on  a  blunt  cone  in  laminar  flow 

2.2  ■  20  percent  station  of  the  initial  length  of  the  cone  in  laminar  flow 

2.3  -  40  percent  station  of  the  initial  length  of  the  cone  in  laminar  flow 

2.4  ■  60  percent  station  of  the  initial  length  of  the  cone  in  laminar  flow 
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2.5  =  75  percent  station  of  the  initial  length  of  the  cone  in  laminar  flow 

2.6  -  90  percent  station  of  the  initial  length  of  the  cone  in  laminar  flow 

2.7  =  Maximum  diameter  point  on  a  blunt  cone  in  laminar  flow 

2.8  =  Maximum  diameter  point  of  a  sharp  cone  for  the  given  cone  half-angle 

in  laminar  flow 

3 , j  where  J  =  1,2,..., 8  =  Turbulent  flow  for  the  same  body  stations 
used  in  laminar  flow 

4,1  =  sonic  point  in  turbulent  flow 

The  total  mass  loss  rate,  W ,  is  calculated  in  the  subroutine  T0MAL0  and  is  the 
integrated  value  of  m(i,j)  over  the  surface  of  the  body. 

In  the  laminar  and  turbulent  flow  regimes,  the  following  notation  will  be  used 
to  identify  the  various  skin  friction  drag  terms.  The  skin  friction  drag  coef- 


f icient 

will  be  denoted  by  1 

CDf  ( i ,  j,  k)  where: 

i  - 

1:  sharp  cone 

i  = 

2:  blunt  cone 

j  = 

1:  turbulent  flow 

j  « 

2:  laminar  flow 

k  - 

1:  ablation 

k  . 

2:  no  ablation 

Sharp  cone  inviscid  pressure  drag  coefficient  for  the  given  cone  angle: 

%  ■  /'  , 

PLEO  1/2  PmVj 

The  forebody  pressure  drag  coefficient,  as  calculated  below,  is  used  in  the  strong 
Interaction  flow  regime,  laminar  continuum  regime,  and  turbulent  continuum  regime. 
For  sharp  and  blunt  cone  configurations,  the  zero  angle  of  attack  forebody  pressure 
drag  coefficient  is: 
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0C  has  the  units  of  degrees  in  the  above  summations. 

The  corrections  for  angle  of  attack  effects  on  forebody  pressure  drag  are 
presented  below. 


[Logio  (0C)1 '  l  Lo«10  <  I «  I)  1’ 

If  K'  <  0.0  ,  then  set  K'  -  0.0  . 
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A  (421  4  «  f  4  j  t  12k)  (0C)‘  ( j  a  |  )  i  U)k 
for  1  0  |  <4° 

for  4°  <  |  a  i  <  40° 


If  Cn  /Crj  s  1.0,  then  set  CD  /CD  =  1.0, 

no  a  o 

Thus,  the  forebody  pressure  drag  coefficient  for  zero  angle  of  attack  is: 
cDp  =  cDp 

O 

and  for  angle  of  attack  cases: 

CI>P  =  CDpo  ,CDa/CD0' 

NOTE:  If  a  rotational  three-degree-of-freedom  trajectory  is  being 

calculated  and  if  the  period  of  oscillation  of  the  vehicle,  t' , 
is  between  the  upper  and  lower  cycle  time  limits  that  are  input, 
the  only  angle  of  attack  correction  made  is  on  the  forebody 
pressure  drag.  The  correction  is: 

cDp  “  cDp  <1+«) 
r  o 

where  a  is  the  integrated  angle  of  attack  calculated  in  the  angle  of  attack 
section  (SUBROUTINE  ROTATE). 

Necessary  for  the  induced  drag  calculations  in  laminar  and  turbulent  continuum 
flow  are  the  induced  pressure  gradient  parameter: 

Fj(K)  »  0.9  -  0.119  sin  Qc  +  0.0108  sin  0C)2 

and  the  wall  temperature  parameter: 


0.968  Tw (2>8> 


To  determine  in  which  flow  or  fairing  regime  the  drag  calculations  are  made,  a 
check  is  made  on  the  values  of  \  and  xj.  The  endpoints  of  the  fairing  regime 
between  free  molecular  transition  flow  and  strong  interaction  flow  are  defined 


by  vi  and  xi|ow*  The  values  generally  used  are  xiu  ■  0.4  and  Xl1()w  ■  0*2  but 

other  values  may  be  used  in  special  cases  to  provide  a  smooth  drag  history. 
Likewise,  the  endpoints  of  the  fairing  regime  between  strong  interaction  drag 


and  continuum  drag  are  defined  by  \'Up  and  Xlow  1116  values  generally  used  are 
Xup  *  6.0  and  xiow  ”  4*0  but,  again,  other  values  may  be  used. 
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1.  If  VI  >  XI  ,  only  the  free-molecular  drag  Is  calculated. 

2  If  XI  <  XI P  and  XI  1  Xl,„,  the  free-molecular  drag  and  the  strong 

interact lon^drag  are  both  calculated  and  the  fairing  technique  for 
values  between  those  two  regimes  is  utilized. 

3.  If  x,  <  xi,ow  and  X  >  Xup*  onlV  stron!;  lnteractlon  drag  iS  calculated‘ 

—  .  “  v  cfmna  interaction  draft  and  the  laminar 

4'  are  both'calculated^and  the  fairing  technique  for  values 

between  those  two  regimes  is  utilized. 

5.  If  x  <  Xlow  and  z  2  ZTR>  the  laminar  continuum  drag  is  calculated. 

Tf  7  c  7  and  Z  >  (ZTB  -  20,547.8  )  both  continuum  laminar  drag  and 
6-  "„tinuuJRturh„lent(dTrRag  are  calculated  and  the  fairing  technique 
for  values  between  those  two  regimes  is  utilized. 

7  If  Z  <  (ZTR  -  20,547.8  )  the  continuum  turbulent  drag  is  calculated. 

The  free-molecular  transition  drag  regime  is  described  by  the  following 
system: 

(Transition  S,»n»)  -  P(.)  <F.«r  Mnin.nl.  Flo.  Sy.r.nO  *  0  -  P(.))(C«tan.  *— > 

.here  P  (c)  1.  the  probability  that  a  molecule  will  collide  with  the  surface 
before  colliding  with  another  molecule. 

A  curve  fit  of  P  (O  as  a  function  of  Knudsen  number  yields  the  following 
relationship: 

2.11053  x  IQ-7  J  1305  r. 

Xw  ”  12.25  p^M*,  f  yTM 


D  -  2  Rn 


zu 

E 

i  -  0 


B(i)  Log10 


G)'] 


P(c) «= 


for  —  >  0.04 

D  " 


0.506  -  0.147  Log10 


/ 0.04  D\ 

\K  / 


for  —  <  0.04 
D 
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For  a  sharp  cone  the  above  system  reduces  to: 


C°c  PWCDr„",-PWCD, 


The  free-molecule  drag  coefficient  is: 


sin  fl  y/W  2M, 


2  \  T, 


2  sin2  0„ 


sin  9.  \  ~n  /  T, 


+  1  + 


I  1  +  erf  (Mm  sin  9  )  ] 


where  fx  is  an  axial  accommodation  coefficient  taken  to  be  0.9.  The  continuum 

drag  will  be  taken  as  the  Newtonian  value  Cn  «=  2  sin2  9,  . 

dN  c 

For  the  blunt  cone,  we  perform  the  analysis  similar  to  a  free-molecule  or  a 
Newtonian  analysis.  The  drag  on  a  blunt  cone 


/  RN  \ 2  2  f  /  rN  \  2 

M*B") 


here  CD  is  the  value  of  the  drag  obtained  for  a  sharp  cone  as  shown  in  the  above 
equation.  Ine  quantity  CDs  is  the  free  molecule  drag  coefficient  on  a  spherical 
nose  and  a  curve  fit  of  experimental  data  yields 


Cfj  *  Cr\ 

°S  DFMg 


A  (174  +  n)  Log  |Q 


r'-ld 


where  the  free-molecular  drag  for  a  sphere  is 


/  -  .  /  .1  - M  2  . 

2*1 - \  -  / — 1— \  erf  (M  )  +  - -  /—  +  —\ 

\*J  j  \mj)  ^  \hJ  m-/ 

\  1  £ 

J  \  T-  / 


Since  the  sphere  drag  is  rather  a  small  portion  of  the  total  blunt  cone  drag  we 
will  approximate  the  above  result  for  hypersonic  speeds 


Cn  -  2.0 
UFMS 
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Thus,  for  a  sharp  cone  in  free-molecular  transition  flow 

CD  cn  C'n 
u  UFM  DC 


and  for  a  blunt  cone  in  free-molecular  transition  flow 


CD  CDcw  =  CD. 


The  drag  in  the  strong  interaction  regime  is  composed  of  interaction  effects  and 
pressure  drag.  The  interaction  drag  is  defined  as  follows: 


/ 


1  1  I 

ZEE 

i  =  0  j  =  0  k  =  0 


A  (200  +  i  +  2  j  +  4  k) 


T„,  l  + 


for  0C  <  15' 


1  1  1 


J  “  “ 

r  i  -  o  j . 


A  (384  +  i  +  2  j  +  8k) 


0  k  -  0 


(l  +  - — -  M  2 

i  2  00 


forOc  >15° 


The  total  drag  in  the  strong  interaction  regime  for  sharp  and  blunt  cones, 
including  angle  of  attack  effects,  is: 


cos  (a  ')  +  Cr 


cn  *  Cr\  »  (Cn  “  Cn  )  1  —  !  — — 

D  Ds  dst  dpleo  V>„ 


The  fairing  between  the  free  molecular  and  strong  interaction  regimes  is  a 
linear  weighting  based  on  the  rarefaction  parameter  value  between  the  drag 
coefficients  of  both  regimes  evaluated  at  the  point  of  interest. 


CD  *  CDt 


~  *»low 
*lup  "  Xllow 


+  Cn  1 


“  Xllow 
Xlup  ■  Xllow 
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For  laminar  continuum  flow,  the  drag  is  composed  of  forebody  pressure  drag,  base 
pressure  drag,  skin  friction  drag,  induced  pressure  drag,  pressure  gradient  in¬ 
duced  skin  friction  drag,  and  transverse-curvature  induced  skin  friction  drag. 

The  forebody  pressure  drag  is  given  in  the  preliminary  calculations  section.  For 
the  base  pressure  drag,  a  curve  fit  of  test  data  gives 


The  sharp  cone  laminar  skin  friction  drag  for  the  no  mass  addition  case  was 
derived  from  the  Blasius  flat  plate  incompressible  solution  modified  for  conical 
flow  and  compressibility.  The  solution  is 


—  =  0.5  +  0.5 
h 


Cp  Tw  (2, 8) 


CD  T. 


+  0.0935  (y-  l)Me 


CDf  (1,2,  2) 


The  mass  addition  correction  to  skin  friction  drag  and  the  boundary  layer 
displacement  effects  in  the  induced  drag  components  are  dependent  upon  the 
laminar  no  blowing  skin  friction  coefficient,  which  is: 
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‘o,  1.53 


1.15  P“lU" 

—  tanfL  - —  Cr>  (1, 2, 2) 


2  uf 


The  blowing  skin  friction  drag  is  then: 

CDf  (1.  2.  2) 

CD  (1,2,1)  =  - 

f  ,  2  m  (2, 8) 

1  +  - 

PCUe  Cf 


The  effect  of  bluntness  on  the  skin  friction  drag,  for  both  blowing  and  no 
blowing  cases,  is  as  follows: 


4  2 

EL 

i  -  0  j  =  0 


A (457  +  i  +  5 j)  [Log10(Re  )]*  — 


L  \  R 


for  -  <  0.2 


CD  (2,  2,  k)  =  CD  (1,  2,  k) 
f  f 


4  2 


i  =  0  j  »  0 


A  (472  +  i  +  5  j)  l  Logj0  (Re  )]*  — 

L  \  rB 


rN  \' 


for  -  >  0.2 

RB  ‘ 


where  the  k  in  the  skin  friction  drag  coefficient  subscript  can  be  set  for 
either  mass  transfer  or  no  mass  transfer. 

If  mass  transfer  and  aerodynamic  heating  effects  are  being  considered  in  the 
calculation  of  the  induced  drag  components,  the  wall  enthalpy  is  set  equal  to 
the  wall  enthalpy  calculated  in  subroutine  EVIL: 

«w  -  ^ 

If  a  nonablating  case  is  being  considered,  the  wall  enthalpy  is: 


Hw  -  0.24  Tw  (2, 8) 

The  recovery  enthalpy  is: 


Hr  -  0.9  H, 
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lor  all  three  components  of  induced  drag  in  laminar  flow,  the  change  in  effective 
cone  angle  due  to  displacement  effects  is: 


»eff-0c 


Cfo  i  2  m  (2,8) _ 1 _ 

—  ^Cf  +  j  + 


1  3 


P e  ue  H. 


•  /  Hw\  ’ 

A  (283  +  i  +  2  j)  (M.2)1  — 


i  =  0  j  =  0 


m  (2,8) 


Hr/  \ 


The  induced  pressure  drag  coefficient  is  therefore: 


CD  =  1.33  \[T  1  - 


cos  0C  CD 


(fleff  “  9C)  Fi  (K)  Me 

The  pressure  induced  skin  friction  drag  is: 


Me  yT 


CD  (1,2,2)  Fi  (K) 


Tw  (2, 8) 

(0eK  -  ec)  -  0,823  +  0.524  — - -  +  0.483  Me' 

L  ie 

The  transverse  curvature  induced  drag  is: 


\c  Me2  tan  dc  V1T 


1  -  (  —  cos  dc  I  Cjj  (1,  2,  2) 


Tw  (2, 8) 


(0eff  -  ec)  0.517  +  0.913 


+  0.0484  Me2 


The  total  induced  drag  for  laminar  flow  is  cDilam  =  cDIp  +  CDIf  +  CDj  * 

Thus,  the  total  drag  coefficient  for  laminar  continuum  flow  is: 

CD  '  C°L»M  ‘  %  +  CdB  +  CD1  2'  k>  *  XaM 

where  the  appropriate  skin  friction  drag  is  used  depending  upon  the  nose  blunt¬ 
ness  and  ablation  rate. 
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The  fairing  between  the  strong  interaction  regime  and  the  laminar  continuum 
regime  is  a  linear  weighting,  based  on  the  viscous  interaction  parameter  value, 
of  the  drag  coefficients  of  the  two  regimes  evaluated  at  the  point  of  interest. 


D  * 


X  ~  X'loi 
X'up  ~  X'loi 


♦  C 


dlam 


For  turbulent  continuum  flow,  the  drag  coefficient  is  composed  of  the  forebody 
pressure  drag,  base  pressure  drag,  skin  friction  drag,  and  induced  pressure  drag. 
The  forebody  pressure  drag,  along  with  angle  of  attack  corrections,  is  given  in 
the  preliminary  calculations  section.  The  base  drag  equation  for  turbulent  flow 
is  the  same  as  the  equation  for  laminar  flow  and  the  appropriate  equations  can 
be  found  in  the  laminar  flow  section.  The  sharp  cone  turbulent  flow  skin  friction 
drag  for  the  no  mass  addition  case  was  derived  from  the  Shultz-Grunow  flat  plate 
incompressible  solution,  modified  for  compressibility  and  conical  flow.  The 
solution  is  as  follows : 


0.5  +  0.5 


CP  Tw  (3, 8) 


Cp  Te 


+  0.099  (y  -  1) 


CpeTe 


3.5964  h* 


for  h*  <  1110  Btu/Ibm 


3  2 


A  (162  +  i  4j) 


i  =  0  j  =  0 


,hv  &y 


for  h*  >  lllOBtu/lbm 


T* 

1  /  Pe  \" 

U 

1  8  L°®10  (2II6  / 

ZZ  =  2.5  +  0.1  tan  h  (  —  -  7 
\  500 


0.4 tan  h  (  — — - 7)  +tanh[— -  5.8 ) 

\  1000  /  \2500  / 
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(2.27  x  10“8)  (32.2)  [  T*  1 1  5 
T*  +  198.6 


fori,  >  2.0  feet 


The  mass  addition  correction  to  the  skin  friction  drag  and  the  boundary  layer 
displacement  effects  in  the  induced  pressure  drag  is  dependent  upon  the 
turbulent  no  blowing  skin  friction  coefficient,  which  is: 
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In  calculating  the  induced  pressure  drag  if  mass  transfer  and  aerodynamic  heating 
effects  are  being  accounted  for,  the  wall  enthalpy  is  set  equal  to  the  wall  en¬ 
thalpy  calculated  in  the  subroutine  EVIL: 

Hw  -  Hw 

If  a  nonablating  case  is  being  considered,  then  the  wall  enthalpy  is: 

Hw  -  0.24  Tw(3,8) 

The  recovery  enthalpy  is: 
hr  *  0.9  Ha 

The  change  in  effective  cone  angle  due  to  displacement  effects  is: 
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f°T  [  2  m  (3, 1 

f>df  -  °c  3.h  '  pe  UeC 

IL 

H.  I 

0.M7  -  *  0.33  *  0.68 

»R  \ 

».\  1 

*  0.106  -  1  M.z  .  +  - 

»R  /  J  3 


PeUeCL  1.2  m  (3.8) 

°T  1  *  - 


*  0.106  —  |  M.2  I  + 


1.6  m  (3,8) 
3-6  Pc  Uc 


The  Induced  pressure  drag  then  becomes: 


cD|p  Ml  V? 


'  £  ”*■)  _ 


I’c  v’<  CL 


0.53  *  0.68  —I  Mc  4  ^0.083 


CD  «?eff  -0c)Fi(K)Mc 
/  pLHO 


and  the  total  induced  drag  for  turbulent  flow  is  cI>iTUR|,“  CDjp  •  T'lus  the  total 
drag  coefficient  for  turbulent  continuum  flow  is: 


C°  C°TURB 


CDp  +  CDh  *  l.k)  4  CD 


where  the  appropriate  skin  friction  drag  is  used,  depending  upon  the  bluntness 
ratio  and  ablation  rate. 

For  the  fairing  region  between  the  laminar  and  turbulent  regimes,  only  the  skin 
friction  drag  and  total  induced  drag  are  modified  because  the  forebody  pressure 
drag  and  base  drag  equations  are  the  same  for  laminar  and  turbulent  flow.  The 
fairing  process  is  as  follows: 

=  . 1L  Pv 

,m  j|_  M  )J  '  L  '• 


(La)  q  (La)' 


The  faired  induced  drag  is: 


<•»,  PCD,  *(1-P)CD 
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and  the  faired  skin  friction  drag  is: 

CDf  **  P  Cpj  (*•  2,  k)  ♦  (1  -  P)cDf(i-  llk) 

Thus,  the  total  drag  coefficient  in  the  boundary  layer  transition  regime  is: 

CD  ”  CDp  +  CDB  +  CDf  +  CDj 
3,1.6  Equations  of  Motion 

Plight  tralectory  points  and  the  associated  flight  parameters  are  computed  by 
ai  iterate  process  in  subroutines  DEREQ,  TEQUAT,  ADM4RK,  and  part  of  ROTATE 
This  involves  the  computation  of  the  altitude  derivatives  of  the  thrusting  and 
tralectory  parameters  at  a  given  trajectory  point  and  the  integration  of  the 
derivativesPover  an  altitude  increment  to  determine  the  parameter  values  at  the 
nex^ tra jectory^oint ,  The  altitude  increment  over  which  the  integration  is 
u7,pL,nt  upon  the  M*nltudo  of  the  p.ranutor  derivative,  and  1, 
controlled  by  the  predictor-corrector  in  subroutine  ADM4RK. 

The  comoutation  of  the  altitude  derivatives  of  the  thrusting  and  trajectory 
par.™«»  uTnl  in  ,ubroutl„e,  DEREQ ,  TEQUAT,  and  ROTATE.  The  proee  are  used 
is  a  compute  the  time  derivatives  of  the  parameters  and  b.  change  to  altitude 

dependent  derivatives  by  means  of  an  appropriate  a™*ck 

necessary  input  quantities  are  the  drag  coefficient,  CD,  total  angle  of  attack., 
instantaneous  mass  of  the  vehicle,  m,  normal  force  coef ficient,  CN  ,  pitching 
“oSn'i“n».cM.CH,,  tho  Euier  anglu.,(«„  ,  4.  #;>.  th.  conpon.nf.  of 
the  rotation  rate  (p  ,  q  ,  O  q,  the  moment  of  inertia  about  the  body  axis,  I  (in 
tabular  form  as  a  function  of  altitude),  the  moment  of  inertia  along  the  body 
axis,  I  (also  in  tabular  form  as  a  function  of  altitude),  and  the  flow  field, 
thrusting,  and  geometry  parameters.  The  drag  coefficient  and  angle  of  attack 
histories  may  either  be  calculated  through  the  appropriate  subroutines  or  be 
input  as  a  function  of  altitude  in  tabular  form.  Any  appendage  on  the  vehicle 
may  be  accounted  for  by  inputing  an  appendage  drag,  C^,  table  as  a  function  o 

altitude.  This  will  be  added  to  the  vehicle  drag  to  give  the  total 
history  which  is  used  in  the  trajectory  calculations.  In  order  to  th® 

input  added  drag  tables  to  be  used  repeatedly  for  different  sired  veh*^c®’ 
input  added  drag  is  based  on  a  constant  input  reference  area,  Awkef  ,  nnd 
lsPscaled  to  the  proper  vehicle  reference  area.  The  total  vehicle  drag  is  then. 


ldt  ' 


(S?) 


Changing  the  thrust  vector  from  body  coordinates  to  trajectory  coordinates 
gives  (lbf): 

Tu  -  T„Y  (cos  0a  cos  dr)  +  T|i  l  co. d-  »•"  0a  sin  d 


-  sin  i/r  cos  (^1  4  lco»  ^  c0*  ^ 


4  sin  0  sin 
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tH  T„  _  (cos  0tl  sin  <Ji )  t  T|1  l  sin  ifr  sin  0Q  sin  4>  -  cos  t],  cos  <f>] 

Y  *jp  X  * 

+  [  sin  i //  sin  0a  cos  <£  -  cos  i //  sin  ! 

tHz  “  T„x  sin  0a  ♦  THy  cos  sin  *  4  THy  cos  0a  cos  d> 

The  velocity  vector  in  trajectory  coordinates  is  given  by  (ft/sec) 

Xg  -  cos  0Q  COS  <!' 

Vp  »  Vm  [  cos  i//  sin  6a  sin  <*>  -  sin  0  cos  1 
Z„  «  |  cos  lit  sin  0a  cos  <A  t  sin  1 1>  sin  <£| 


The  normal  and  axial  force  coefficients  are: 


Fn  «  CN 
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ARef  Cx  sin  Oa  - 


'^B2  *  ^B2 


cos  &a  sin  <{> 


V  2  .  y  2 
YB  +  /jB 


COS  0a  COS  0 


The  altitude  derivatives  of  the  thrusting  and  trajectory  parameters  are  calcula¬ 
ted  by  multiplying  the  thrusting  and  trajectory  time  derivatives  by  the  derivative 
of  time  with  respect  to  altitude  (dt/dz).  The  derivatives  calculated  in  TEQUAT 


Velocity  Derivative: 


dz  U  sir 


™*i  P-uJ  cDARef 


-  32.21852 


Re  +  Z 


where  R*  is  the  radius  of  the  earth. 


Flight  Path  Angle  Derivative: 
dyf  ,  ^  cosy, 

dz  sin  y,  J  Um 


Fz  +  Thzt  I 


U«2  /  Re  \2" 

-  -  32.21852  - ) 

Re  +  Z  \Re  +  Z/ 


Range  Derivative: 


lxR  r  /  Re  \  "1  1 

-  •  cos  (dr  )  cos  yt  [  -  I  IL - 

dz  Yt  \  Re  +  Z  /  w  U,  sin  y, 


Cross  Range  Derivative: 


dYR  1  T  .  /  Re  \  I 

-  „  -  (in  k  j  cos  yt  I -  I  I 

<1.  V„.inr,  [  w->  "  U,.z/  J 

Thrust  Misalignment  Angle  Derivative  In  The  Yaw  Direction: 


U^siny,  LmU^cosy,  J 
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Thrusting  Weight  Derivative: 


>'f  V  *SP 


The  derivatives  calculated  in  ROTATE  are: 

Euler  Angle  Derivative  In  The  Yaw  Direction: 


d*  sin  Yf 


R  cos  0  »  Q  sin  0 


Euler  Angle  Derivative  In  The  Pitch  Direction: 


dz  U  sin 


(Q  cos  0  -  R  sin  0) 


Euler  Angle  Derivative  In  The  Roll  Direction: 


d0  1 

dz  sin  y{ 


P  4  U,K  sin  Y(  ~  'in 


Pitching  Component  of  Rotation  Derivative: 

dQ  1  \  P".  u-2  Rfl  AKcf 

-  m  - . — -  . - — -  Cu  (sin  i 

dz  Uw  sin  yf  I  1 


RB  CM„ 


Cy  (sin  6a  cos  0  cos  0  sin  0  sin  0) 


(^) 


Ablation  Weight  Derivative: 


dz  U_,  sin 


Yawing  Component  of  Rotation  Derivative: 
dR  1  U~2  RB  ARe(  T 

—  «  . . -  . . . . . -  Cj^|  (sin  cos  <5  -  sin  0a  cos  i{/  sin  $) 

dj!  Vm  Yl  |  1 


rb  cm„ 


*zl 

j  '  \  I  | 
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Rolling  Component  of  Rotation  Derivative. 


dz  sin  y(  V  I**  J 

,  int-PoratinE  the  above  derivatives 

The  next  trajectory  point  is  det“'"i^  alte^s  the  range  (AZ)  over 

by  means  of  the  predictor-corrector  method  which  al  vcls  are  maintained. 

which  the  integration  is  technique,  it  will  not  be 

Since  the  predictor-corrector  method  rical  description  of  the  program 

explained  here,  but  can  be  found  in  t  velocity  at  the  next  trajectory 

under  the  ADM4RK  subroutine.  Then,  simply,  the  y 

point  is: 


.Z-  AZ 


1*.  raining  toi-tiv-  „  t-JJ..- J  l^he 

L  vlght  lo..«  «.-  thrusting  .nd 
ablation  over  the  integration  altitude  interval. 


That  is: 


/Z  -  AZ 

(S)* 

£. 

These,  then,  «  —  to  cl,...-;.  .-—  -***  **  ''‘hl'U'  “ 

stated  In  the  preliminary  calculation  section 


W  -  W0  AW  -  AW-™ 
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3.2  WAKE  CALCULATIONS 


The  decoy  design  program  requires  calculation  of  the  radar  cross-section  (RCS) 
and  the  wake  length  as  a  function  of  the  decoy  configuration.  This  involves 
wake  configuration  calculations  to  define  the  wake  flow  field  properties,  and 
the  RCS  and  wake  length  calculations.  The  input  quantities  necessary  for  the 
wake  configuration  calculations  are  the  flow  field  and  body  geometry  properties, 
transition  altitude,  ZTR,  drag  coefficient,  CD,  and  the  flight  path  angle  ye  . 

Also  needed,  and  available  through  input  tables,  are  the  heat  shield  conductivity, 
(Btu/f t°-R-Hr) ,  heat  shield  specific  heat,  Cpw  (Btu/f  t°-R-llr)  ,  chemical  enthalpy 

of  the  heat  shield  AhCHE!^ft2/sec2) ,  heat  shield  density ,  pw  (lb/ft3) ,  heat  shield 

thickness,  8  .  (in),  ablation  temperature  of  the  heat  shield, TA(U  (°K),  sea  level 
density,  (lb/f t3) ,  scale  height,  , (kft) ,  and  scaling  constants,  bn  ,  b21  ,  b22 , 

b23>  b2i • 

3.2.1  Preliminary  Wake  Calculations 

The  calculations  are  as  follows: 

Altitudes  are  used  in  terms  of  thousands  of  feet. 

ZTR 


Z 


Z 

1000. 


The  total  drag  area  is  given  by: 

cDt  A  .CD(rrRB2) 

The  viscous  drag  area  by: 

CDV  A  -  <CD  -  <CDp  +  cDb  *  CD,  >  1  *rB2 
The  velocities  have  units  of  thousands  of  feet  per  second. 

U  *  ■■■- . 

"  1000. 

Ue 

n  .  - 

e  1000. 


The  mass  ablation  rate  is: 


m*ABL 


-  W 
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The  mass  swallowed  by  the  boundary  layer  is: 


v  0.25 


0.5 


M. 


,0.5 


2010  \7 


(*R> 


1.5 


U  M. 


0.5 


^OO  -Xi  *00 


pl 


sin 


mBLS 


for  1.98  Rg  >  2  Rj^ 


C76  P<»  u~  "■  rB 


2000  px  UM  RB 
M„ 


-77 


for  2  Rfj  >  1.98  Rg 

The  length  of  the  conical  frustum  along  the  cone  is: 

0.0 

for  2  Rpj  _>  1 .98  Rg 

Rg  -  Rn  tan  6C 
for  1.98  Rg  >  2  Rn 
Base  diameter  is: 

Dg  =  2  Rg 

3.2.2  Flow  Field  Calculations 

The  wake  flow  field  calculations  are: 

Free-Stream  Reynolds  Number  Based  on  the  Body  Diameter 


Re 


UD 


10 3  P„  Dg 


Kw  sin  ye 
Pw  CPw  ®2wH 
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Conical  Wall  Temperature 


278.  f 


C6o  (Ky  Z) 


'59 


i  z  >  ZBLT 


*wc 


Spherical  Wall  Temperature 


;  z  <  z 


BLT 


T  W8P  ~  278<  4  C60<KvZ)C59 

\  Kw  DB1/2  exP  (Z/4^z^ 

(  ^3/2  swH  \ 

Total  Mass  Rate  in  the  Boundary  Layer 

.  *  ■  *  •  * 
m  BL  “  m  BLS  +  m  ABL 

Wall  Temperature 

l  TABL  •  (TABL  i  TWC  and  2RN  ^  °-"DB>  ot(TABL  <  Twsp  and  2RN  >  0.99  D0) 

Tw  =  Twc  *  2RN  <  °-"  db  and  Twc  <  tabl 

f  TVSP  ’  2  rn  -  DB  and  TWSP  <  TABL 

Wall  Enthalpy 

=  1.087  (104)  T„, 


Drag  Area  of  the  Second  Entropy  Layer 


ct)vS2  A  =  CdtA  “  C°vA 
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Oq  «=  max  ( Oc  ,  0C  ) 

Soc  c 

(9-  =  min  (0c  >  90-) 

Soc  ^ 

Mass  Flow  Rate  in  Second  Entropy  Layer 
m^s 2  =  2000.0  p„  UK  RN2  cot2  0 ^ 

Velocity  Along  the  Cone  in  the  Second  Entropy  Layer 


500  px  U J  CDis2  A 


U2c  =  max  (U2c  ,  1.0) 


A,  =  0.286  +  1.029  M 


Maximum  Shock  Angle 


0<-2  =  57.3  sin' 
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TABLE  II 1-2 
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Number  of  Electrons  Produced  by  the  Nose  Cap 
(30.48)  ^ 


Cone  Pressure  with  Respect  to  Second  Entropy  Layer 
(7Mj  sin ^  fU  -  1)  P 


Cone  Enthalpy  with  Respect  to  Second  Entropy  Layer 

K  P2c  <p2c  +  6  p~  > 

'2C  p.  (  «  p2e  ♦  K  ) 


Cone  Density  with  Respect  to  Entropy  Layer 


P2c  *  6  K  1 


*p~_ 


Definition  of  the  Flow  Properties  at  the  Cone 

If  mBLS  ^  "‘jsj 
Pc  m  p  f 


hc  -  he 


pc  *  (p«/p~)  p- 


Uf  «  ue 


If  m  bxs  <  m  ls2 
pc  •  pjc 


hc  -  h2c 


pe  *  P2e 


Uc  *  U2< 
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Chemical  Length 

I  C 
P2  c 


117 


Sc  mlS2 


„  17 

pt  ~c 


(m  BLS  “  m  iS2  >  .  , 


lJ2e  "'*BLS 


We  m  BLS 


:  m  BLS  mXS2 


XHEM 


C 
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f*2c  Sc  m  BLS 


m 


BLS 


<  m\ 


S2 


U,„  m 


LS2 


Number  of  Electrons  at  Shoulder  of  Nobc  Cap 

Ci  i2  NRN  m*N  (I  -  Cj2j  +  (1  *  Nrn  Schem)  Cm  e*p(-bJ5  SCHEM)  ] 


N, 


SRN  ” 


N  +  C123  SCHEM  NRN  TW 


124 


wHI 


|  Cc  hc  } 

h  < -  ■  -rr-  •  mrat 


I  pl 


RT„ 


RAT  j  (TABLES  111-3  THRir-7) 


n«BL  "  nfd*H 


118 


i  -  eip  (-  b22  Schem  ) 


"1 

1 

/  Uc-22\' 

b2l  *  b2)  1 

(  »  jj 

Number  of  Electrons  Leaving  Boundary  Layer 


(30.48) J  neBL  B  BL 


^BL 


Number  of  Electrons  Entering  Wake  Neck 
*  NBL  *  NSRN 

Cone  Hach  Number  with  Respect  to  Second  Entropy  Layer 


M2c  -  kmi  ( M2c ,  l ) 


Second  Entropy  Layer  Hach  Number 
M2  «  I,  m2e  ,  \  )  (Table  I1I-8) 
Wake  Mach  Number 

%  -  »,  (M#.  0C>  (Table  111-8) 


# 
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TABLE  U»-3 

9#  AS  A  FUNCTION  OF  NORMALIZED  DENSITY  AND  ENTHALPY 


TABLE  II 1—4 


ELECTRON  DENSITY  AS  A  FUNCTION  OF 
NORMALIZED  ENTHALPY  AND  AIR  DENSITY 
FOR  1000  ppm  SODIUM  SEED,  mrat  -  0.0 


n 


t 


p/pa  >  [clcrtrons/cc] 


h  y 

^'PolATMl 

Creator 

than 

0.1 

0.1 

0.03 

0.01 

Leaa 

than 

U.01 

Less  than  IS 

0.0 

0.0 

0.0 

0.0 

0.0 

IS 

0.70x10* 

0.  50x10* 

0. 32xl02 

0.95xl02 

0. 90xl02 

20 

0.38xl03 

0.38xl03 

0.12xl04 

0. 30x1 04 

O.lOxlO4 

40 

0.32xl010 

0. 32xl010 

O.ilxlO10 

0.20xl09 

0.20xl09 

60 

0.25xl012 

0. 25x1 032 

0.78X1011 

0.05X1011 

O.OSxlO11 

80 

0.32xl013 

0.32xl013 

O.lOxlO13 

0.14xl012 

0. 14xl012 

100 

0.17xl014 

0.17xl014 

0.6Sxl013 

0.20xl013 

0.20xl013 

120 

0.84xl014 

0.84xl014 

0.31xl014 

0.78xl013 

0. 78xl013 

140 

0.24xl015 

0.24xl035 

0.87x10** 

0.22xl014 

0.22xl014 

160 

0.44xl015 

0.44xl035 

0. 16xlO*3 

0.38xl014 

0.18xl0*4 

Craatar  than  160 

0.44xl015 

0.44xl033 

0.16xl015 

0. 38xl014 

O.lSxlO14 
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TABLE  lll-S 

ELECTRON  DENSITY  AS  A  FUNCTION  OF 
NORMALIZED  ENTHALPY  AND  AIR  DENSITY 
FOR  1000  ppm  SODIUM  SEED,  mrat  ■  0.0) 

ne  -  nf  j  A-  .  p/p0l  Ulcctfon./cc) 
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TABLE  III-6 


ELECTRON  DENSITY  AS  A  FUNCTION  OF 
NORMALIZED  ENTHALPY  AND  AIR  DENSITY 
FOR  1000  ppm  SODIUM  SEED,  ^rat  -  01 


TABLE  111-7 


ELECTRON  DENSITY  AS  A  FUNCTION  OF 
NORMALIZED  ENTHALPY  AND  AIR  DENSITY 
FOR  1000  ppm  SODIUM  SEED,  MRAT  -  1.0 


ne  -  nf  j  ,  p'p0  J  [elcctrons/cc] 
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TABLE  OF  M  VERSUS  bAQ  AND  0Q 


MC 

0 

2 

6 

10 

14 

18 

22 

26 

30 

45 

>45 

1.0 

1.0 

1.13 

1.29 

1.43 

1.57 

1.84 

2.00 

2.13 

2.76 

100 

3.0 

3.0 

3.11 

3.33 

3.58 

3.85 

H 

4.49 

4.88 

5.32 

7.79 

100 

6.0 

6.30 

7.00 

7.84 

8.89 

12.0 

14.6 

18.4 

100 

100 

9.0 

9.64 

11.2 

13.4 

16.7 

21.8 

31.6 

53.0 

100 

100 

100 

12.0 

12.0 

13.1 

16.2 

21.0 

29.8 

51.0 

100 

100 

100 

100 

100 

15.0 

15.0 

16.8 

22.0 

33.4 

58.0 

100 

100 

100 

100 

100 

100 

18.0 

18.0 

20.6 

29.0 

49.2 

100 

100 

100 

100 

100 

100 

100 

21.0 

21.0 

24.6 

37.7 

80.0 

100 

100 

100 

100 

100 

100 

100 

24.0 

24.0 

28.8 

48.6 

100 

100 

100 

100 

100 

100 

100 

100 

27.0 

27.0 

33.4 

63.0 

100 

100 

100 

100 

100 

100 

100 

100 

30.0 

30.0 

38.0 

81.0 

100 

100 

100 

100 

100 

100 

100 

100 

>30.0 

100 

100 

100 

100 

100 

100 

100 

100 

100 

100 

100 

III— 95 


Us  -  U2 


Shoulder  Reynolds  Number 


Ps  Ms  ReUD  C169 


cs  =  c  1 64  +  C165  sin  °c  Re 


in 1/4  0.  Re..  _3/8  M  ~l/2 


Neck  Enthalpy 

hN  =  K  +  CS  [  hs  -  h*  +  0.5  (106)  Us2 


nCHEM 


hn/hs 


2  hs  +  106  Us  2 


Electron  Density  at  the  Neck 
Ps  Ns  hs 


(30.46) 3  m*BL  hN 

/cdta  p.  \1/2 


B  LS  = 


Shoulder  Wake  Momentum  Thickness 


9BLS  >  m*BLS  ^  m£S2 


Cn_ .  A  + 


■"‘2S2  cD2S2a 


’>  m*BLS  <  m*2S2 


Shoulder  Reynolds  Number  Base  on  Wake  Momentum  Thickness 


0%  Res 


hv  1  - 


2  dO6)  Cn0  ps  0*s  Us2  Re0s  (1  +  Hn/Hs) 


*  z  >  ZBLS 


0.0  ;  Z  <  ZBLT 
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Transition  Electron  Density  Scaling  Factor 


b5  =  Cgj  1  +  Co,<  Md  it  ^ 


84  MRAT  +  C86  mRAT  +  (  8 


bll  '  C91  +  c92  MRAT 
Electron  Decay  Rate 


C69  C1 15  +  bll 


'°6  “s2  «BLS 


K«  ■  («,•)' 


7.5  (1028)  pg  Cj 


fell 


4.0  C90  hs  K, 


'N 

e  PS  K'ds 


I  Z  )  Zb 


1.0  ;  Z  <  ZE 


Transition  Electron  Density 

(  /l0~3  °124  -  bt  0  c 


'  I 


C69  b5  Ps  133  gt\  f  (  / 

77^ -  1-0-Mp)"'hN'(cl3’  + 


C136  M  130  +  c131  !hN  I  Cl32)|  + 

1,11  (1°27)  P*  Cl0°  +  "«N  GTCHEM  ' 


AC0N 

^ON 


3 ■  2 • 3  Radar  Cross  Section  and  Length  Calculations 

?'  ln;ut*  for  the  r»d*r  cross  section  (RCS)  and  sake  length  calcula- 

frequency  a  (4  ?'  "ul"  '  <  Men),  see  level  colileion 

xrequency,  v.lCPS),  radar  noise  level  to  which  the  wake  length  is  measured 

“NL  m  )»  altitude,  z  (left),  free-stream  Mach  number,  M„  scale  height  tfv  (kft) 
look  angle  as  a  function  of  altitude,  *  (degrees) ,  el.ct  on  dec.v  r'te  bf  ' 

d'n,lty'  »"  S£S  .‘.-the 


Dw  -  D 


K  ■=  2  it  f  2  x  10^ 
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Surface  Scattering  Correction  for  Base  Diameter 

i),)' 

Dw*  — - - 

2*la2  *  *  h*  l)«> 

Atmospheric  Density  as  a  Function  of  Altitude 

7. 

h. 

/>»  *  * 


Electron  Density  at  which  the  Plasma  Becomes  Ovcrdenae 


’•<  «  ’  fe)  _'  °  •  (tt) 

Transition  Point  Onset  of  Ovcrdcnso  Region 


i  (“  -  777) 

Onset  of  Constant  Multiple  Scattering  Region 

\^)u  fe) ! 1  *>-> ix<* 

*OD  •  \  '  '  ’ 

j  *1  BOO  •  *2  BOO  ‘  X»B 


^  XOD  1  X00  i  °‘ 

/  0.0  ;  XqO  *  0 
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Point  of  Onset  of  Single  Multiple  Scattering  Region 


XSSP 


‘*55 


%S  1  XSS  - 

f  0.0  ;  %  <  0. 

Mallei  of 


3.r  a  10J 
2.  co*  <A 


Point  of  Onset  of  Variable  Multiple  Scattering  Region 


\  XMS  •  XMS  t  XOOP 
■  |  *OOf‘  XMS  <  *ow 


X*«P  -  *ODP 


IV  ak  Radar  Cross  Section 


The  wake  length  Is  baaed  on  a  unimodal  function  of  the  axial  distance  In  the  wake, 
"  (*>  ,  and  an  Input  radar  nolae  level,  os  L  .  The  wake  length  is  the  distance  from 
the  maximum  value  of  e(i)  to  the  point  where  the  function  Is  equal  to  the  nolae 
level.  The  procedure  employed  Is  baaed  on  a  course  marching  procedure  to  define 
the  limit*  for  the  maximum  of  the  function,  followed  by  a  fine  resolution  search 
for  XB>, ,  using  a  Fibonacci  search  technique.  If  the  maximum  of  the  function  Is 
above  the  noise  level,  the  course  starching  procedure  i»  again  initiated  from  the 
peak  and  the  limits  for  the  point  X(ao.  ,  where  the  function  Is  equal  to  the  noise 
level,  are  established.  The  square  of  the  difference  between  the  function  and 
the  noise  level  is  then  minimised  using  the  Fibonacci  technique  to  define  the 
root  more  accurately.  The  wake  length  Is  then  calculated  by: 

*'f  *»oot  “  *«n*i 

The  function,  «(«>  ,  is  defined  by  the  following  equations: 
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150  r 


m 


1  0.0  ;  X  -  X,  0.0  of  if  X  -  X.  > 

\  '  ‘cos* 

/  Vj  Dw*  sin  2  <5  ;  otherwise 


*  X. 


OOP 


Yi  *  smallest  of 


(150  r 
cos  <5 


150  r 


2  LMS  •  X  "  Xt  "  XODP  >  *  XMSP  ♦  -  x 


150  r 


0.0  ,  X  -  X,  <  XODp  or  if  X  -  X,  >  — — -  *  X^ 


Oi  • 

/ 


Yj  Dw*  sin2«J;  otherwise 


150  r 

l  xmsp  ;“^TZ  -  (X"X<  -  xmsp> 


150  t 
cos  6 


-  X,  ;  otherwise 


I  X  -  X,  i  X  <  x,  *  Xssp 

*32" 

/  xssp  ;  x  >  xt  +  x^p 

1 50  r 

If  X-X,  <  x^jp  j*  if  X-X,  >  ~  ^  +  XSSp 


sod 

0.0  and 
<  0.0 
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0.0 


no  t 


if  x  -  x«  £  xmsp  *ml  x  “  xt  7777  *  xssp  ,nd  b:o  *  00 


*3 


fin  ^  <3 


l 

_ 


~)  (■""" 

eCR 7  \ 


bl  V32  \  I 


“20 


b]  b20 


*  2  “TT 


eCR 


bj  b>0 

V,, - Y 

2  41  2 
e  -  c 


J2 


*  b2/— ) 

V 10 «/ 


1.0 


*in‘  d 


(2-b2„)' 


bl  b20 


nor 


If  X  -  X(  >  XMSP  *n<^  x  ~  x.  £  .  ♦  xs5p  *nd  b2o  »  0,0 

co*  *5 


D*‘  2  a  1  /"•*  \  /“bt  Y»l  "bl  M 

ff,  -  »in  ©  —  I  II*  -  e  1 

b>  l  \  Vr/  \  / 

-feK- . ■ 


”  4  bt  VJ2 


(Y}2  -  Yj,)  (tin4*) 


>&)’  - 


1 1.0  4  b2/-i— \  Dw2* 

“  ’Vio*/  *  I 


I 


150  f 


Ul 


\  xsspi  77^  2  (X-  x«  -  xssp> 


1 50  r 


/  X  -  -  -  Xt  ;  otherwise 

eo*  6 


Y42  -  X  -  xt 
If  X  -  X{  <  Xggp 
oA  •  0.0 
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3.3  MISCELLANEOUS  CALCULATIONS 

Four  miscellaneous  operations  are  performed  near  the  end  of  the  nulyols  calcula¬ 
tions.  These  arc:  a.  the  evaluation  of  a  polynomial,  b.  the  calculation  of 
a  vector  of  generalized  differences,  c.  the  calculation  of  the  average  densities 
of  the  decoys,  and  d.  the  calculation  of  a  series  of  terms  comparing  the 
characteristics  of  the  decoys  immediately  before  and  after  the  discontinuous  shape 
change. 

The  polynomial  evaluation  was  originally  intended  for  comparison  of  the  nominal 
free  space  radar  cross  section  of  a  reentry  vehicle,  J[C7*  with  a  curve-fit  of 
the  nominal  radar  cross  section  of  various  decoya.  The  user  supplies  curve  fit 
coefficients  for  the  cross-section  data  applicable  to  the  class  of  electromagnetic 
devices  being  used.  The  difference  between  the  reentry  vehicle  and  decoy  could 
be  calculated  for  use  in  the  penalty  function  equation  as: 


ARCS 


N,  N2  N, 


*C0E(<K- I)  NjN2  *(J-l)N|  ♦ 


I-l  J-l  K- 
RB,  HN, 


This  calculation  was  generalized  to  allow  tha  polynomial  to  be  a  function  of 
any  three  quantities  in  the  0CCUR  vector,  Qj  ,  and  is  now  coded  in  the  following 
form: 
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N, 


'Vcs 


E  EE' 

K  -  I  J  -  I  !  «  1 


OF. 


((1C  -  I)  NjN>  *  ( |  —  I >  Nj  .  I) 


V| 


J-  I  K -  I 

<h  9* 


This  equation  allows  the  radar  calculation  described  above  to  be  perforated 
(  X'(  H  ■  1.0)  and  also  allows  direct  use  of  the  polynomial  as  a  constraint  if 
X(.  is  0.0  and  X(  ^  Is  -1.0. 


An  many  as  twenty  differences  between  pairs  of  quantities  in  the  0CCUR  vector 
can  be  calculated: 


Gj  -  On,  -  Q|.j  l-l.  JO 


where!)  and  1^  are  subscripts  Identifying  the  desired  quantities  in  the  0CCUR 
vector. 


The  overage  density  of  the  Initial  decoy  is  calculated  from: 
r  «  R.^  co* 

h  «  Rj^  ( 1  -  sin  ) 

V,  .  /rh2(RNj  -  h/3)  *("/J)<lA,  -  h)  «B  2  +  Rp ^  r  *  r2  ) 

Dj  -  V,/f, 

If  the  vehicle  undergoes  a  discontinuous  shape  change,  then  the  average  density 
of  the  decoy  immediately  after  the  shape  change  is: 

i  ■  Rj^  co* 
h  -  RNj  (1  -  sin  02  ) 
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V,  .  »h‘<RNj  -  h/H  *(»/J)<La2  -  h)  (RBj  ♦  H:'  *  '2» 


Dj  V2/*2 


The  coopar Isons  of  the  vehicles  before  and  after  discontinuous  shape  change  are 
calculated: 


n,  -  t. 


l)i  ■  1^1  • 

1  1  FINAL  2 


nj  «  *N, 


‘FINAL 


d4  -  rd 


'FINAL 


D«  ■  Ai  ■■ 

’  1  FINAL  2 


°6  *  LA, 


FINAL 


R,  -  f2/», 


FINAL 


FINAL 


f*2'**l 

>.0  if  rNj 


FINAL 


FINAL 


R4  -  *b2'*b, 


FINAL 


(  Kj  !  A  * 

J  2  1  FINAL 

I  0.0  if  A,  -  0.0 

1  FINAL 


R(J  m  L^2  / 


FINAL 


These  quantities  are  available  for  use  as  constraints  (see  Table  4  of  Appendix  II, 
the  User's  Manual). 
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4.0  COMPARISON  OF  DECOY  WITH  REENTRY  VEHICLE 


The  basic  analysis  calculations  described  in  Section  3.0  primarily  involve 
the  determination  of  data  for  a  single  vehicle.  In  this  section  the  techniques 
used  to  compare  a  reentry  vehicle  and  decoy  will  be  discussed  along  with  the 
definition  of  the  corridor  and  effectiveness  integrals.  These  integrals  are 
calculated  for  ubo  in  the  penalty  function  equation  and  in  the  effectiveness 
operations,  respectively. 

4.1  DIFFERENCE  EVALUATIONS 

As  many  ns  nine  performance  functions  are  stored  for  the  reentry  vehicle  and  the 
decoy.  The  differences  between  the  reentry  vehicle  (R/V)  performances  and  decoy 
(D)  performance  at  the  M  altitudes  are: 


U>i  v  Vv,  -  V/*D, 


V*j  R/V:  "  0D; 


*«L.  fL.  '  WL. 

'i  ‘R  Vj  'D, 


-»L, 

2i  ‘D. 


\*L  -  W,  - 

\  '%/Vj  *I>, 


“*>i  ’"'.  v  ■  *% 


'•r,  ■  fRj  “  *R, 


\« 


R.  "  *R. 


’R/Vj  JD, 


These  nine  vectors  of  differences  will  be  referred  to  as  APj  in  the  following 
discussions. 

In  the  influence  calculations  associated  with  M0DE  -  Z  in  the  program,  the 
partial  derivatives  of  the  differences  listed  above  with  respect  to  the  design 
variables,  X,  can  be  evaluated  using  finite  differences.  The  first  derivative 
is  based  on  a  basic  decoy  (1)  and  a  comparison  decoy  (2): 
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The  second  derivative  is  based  on  three  decoys  which  are  used  to  define  two  first 
derivatives.  The  second  derivative  is: 

/<*2\P\  <*'pi  «>i|,  *  <^APj/«?X),  |? 

W'J,  '/»»'*  ~^X2 ,  XX 

4.2  INTEGRALS  OF  SPECIAL  FUNCTIONS 

The  corridor  integrals  are  defined  in  terms  of  the  performance  dif ferences ,  AP  ; 
the  corridors,  T  and  B;  and  an  arbitrary  multiplier  which  can  be  a  function  of 
altitude,  7. ,  if  desired.  The  quantities  AP,  T,  and  B  are  all  functions  of 
altitude.  For  each  of  the  nine  performance  functions  there  is  a  corridor  integral 
of  the  form: 


i  -  1,  M 
i  1.  9 


where  the  square  brackets  mean  that  negative  numbers  are  set  to  zero. 

The  effectiveness  integrals  are  calculated  using  the  performance  differences,  AP  j 
and  the  standard  deviation  of  the  perceived  error  (uncertainty)  in  the  defense 
measurement,  o,  .  Both  of  these  quantities  are  functions  of  altitude.  For  each 
of  the  nine  performance  functions  there  is  an  effectiveness  integral  of  the  form: 
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5 . 0  EFFECTIVENESS  MODEL  OPERATIONS 


A  discussion  of  the  effectiveness  model  operations  is  included  in  the  main  portion 
of  the  ADTECH  IV  Final  Report.  The  effectiveness  integrals,  Ej  ,  are  discussed 
in  Section  4.2  of  this  volume.  For  example: 


2 

i 


The  difference  in  the  means  of  the  logarithm  of  the  likelihood  ratio,  //  ,  is  cal¬ 
culated  using  only  the  terms  which  have  been  specified  by  the  user: 

(Sv  hv  +  SDED  +  S/3E(9  +  SRj  ER2  +  SR2ER2  +  SR3ER3  +  SW11%1+SW2  eW2+  SW2%3)j 


Zo-Zf 


With  the  input  probability  of  false  dismissal,  PpD,  the  upper  limit  of  integra¬ 
tion,  t  ,  is  calculated  from  the  equation: 


I 


dt 


Then  with  tQ  and  /x  calculated  the  probability  of  discrimination  of  the  decoy  is: 


—  oo 
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6.0  CLASSIC  FUNCTIONS  FOR  TESTING  THE  OPTIMIZATION  TECHNIQUES 


Five  functions  or  set  of  functions  are  provided  for  use  in  checking  the  search 
and  optimization  techniques.  The  substitution  of  analytic  functions  in  place  of 
the  complete  trajectory  calculations  significantly  conserves  machine  time  and 
also  provides  positive  control  for  the  checkout  since  the  characteristics  of  the 
functions  are  explicitly  defined. 

The  subscripts  of  the  quantity  X  in  the  following  equations  refer  to  the  locations 
in  the  vector  called  0CCUR  in  the  program.  The  vector  of  coefficients  called  A 
in  the  equations  is  the  same  as  the  vector  of  coefficients  called  A  in  the  drag 
calculations,  thus  it  is  recommended  that  these  calssic  check  cases  and  actual 
trajectory  calculations  not  be  utilized  in  the  same  job  unless  careful  attention 
is  given  to  the  A  coefficient  vector  and  the  0CCUR  vector  for  each  case. 

The  first  function  is  based  on  the  one  used  by  Rosenbrock  (Refs.  III-5  and  III-6). 

,  2  , 

X i oq  =  100  [Xt  —  Xj  ]  +  (1  —  Xjl"  +  aj 


The  second  set  of  functions  provides  two  additional  equations  which  can  be  used 
to  generate  nonlinear  constraints  for  use  with  the  Rosenbrock  function: 


Xioo  =  100  (X2  - 


X/) 


(t  -  X.)' 


X101  =  a2  |X2  -  ^a8  ^Xl  _  a7^  +  “9  -  a7>2  +  ajQ  (Xj  -  a 7)  +  au)| 

X;o2  •  aA  |x2  “  (al3  <Xj  -  *|2>^  f  ai4  (X]  “  ai2^  +  *15  ~  a12>  *  *16^ 


The  third  set  of  functions  is  based  on  the  four-variable  function  with  constraint 
functions  which  was  developed  by  Rosen  (Reference  III-7). 


X200  =  Xt2  +  X22  +  2X32  +  x42  -  5Xj  -  5X2  -  21  X3  t  7  x4  +  ap 


X201  =  -X12  “  X22  _  X32  "  X42  -  Xj  +  X2  -  X}  +  X4  +  8 


x7 


X,2  -2X,2  -X,2  -  2Xx2  +  X,  +  X,  +  10 


202  -  ~  A1  "  /A2  "  A3  “  ^A4  +  A1  +  A4 


x’7q3  =  —  2  X  i 2  —  X,2  —  X,2  —  2X|  +  X7  +  Xa  +  5 


kl  "a2  “  a3  _'a1  +  A2  +  A4 


The  fourth  function  is  based  on  a  two  variable  function  which  was  Fiacco  and 
McCormick  (Reference  III-8). 


X300 


«*1 


■f 


♦  a26 


Ill-Ill 


The  fifth  function  is  a  general  quadratic  function: 

*400  “  31  xr  *  a2  X1  x2  4  “3  x23  +  a4  X1  4  a5  x2  4  a6 

lhis  function  was  used  to  generate  the  optimizer  check  cases  which  were  presented 
in  Section  6.3  of  the  Users  Manual  (App.  II). 
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APPENDIX  III-l 


REPRODUCTION  OF  VARIABLE  METRIC 
METHOD  FOR  MINIMIZATION 


A  report  written  by  William  C.  Davidon 
Argonne  National  Laboratory,  ANL-5990  Rev. 
University  of  Chicago,  November  1959 


VARIABLE  METRIC  METHOD  FOR  MINIMIZATION 


William  C.  Davidon 


This  is  a  m.v  ihod  for  dele  rm.inng  :m:nc  rically  iocal  mininia  of  dif¬ 
ferentiable  functions  of  several  variables.  In  the  process  of  location  ea cr. 
minimum,  a  matrix  which  characterizes  the  behavior  of  the  function  about 
the  minimum  is  determined  For  a  region  ir.  which  the  function  depends 
quadratically  on  the  variables,  no  more  than  N  iterations  are  required, 
where  N  is  the  number  of  variables.  By  suitable  choice  of  starting  values 
and  without  modification  of  the  procedure,  linear  constraints  can  be  imposed 
upon  the  variables. 

IN0T.REPR0DUCTPT" 

1.  INTRODUCTION 

The  solution  to  many  different  types  of  physical  and  mathematical 
problems  can  be  obtained  by  minimizing  a  function  of  a  finite  number  of 
variables.  Among  these  problems  are  least- squares  fitting  of  exper imer.ta - 
data,  determination  of  scattering  amplitudes  and  energy  eigenvalues  by 
variational  methods,  the  solution  of  differential  equations,  etc.  With  the 
use  of  high-speed  digital  computers,  numerical  methods  for  finding  the 
minima  of  functions  have  received  increased  attention.  Some  of  the  pro¬ 
cedures  which  have  been  used  arc  those  of  optimum  gradient, conjugate 
gradients,^)  the  Newton-Raphson  iteration, (3 )  and  one  by  Garwin  and 
Reich.  (4)  In  many  instances,  however,  all  of  these  methods  require  a  large 
number  of  iterations  to  achieve  a  given  accuracy  in  locating  the  minimum 
Also,  for  some  behaviors  of  the  function  being  minimized,  the  procedures 
do  not  converge. 

The  method  presented  in  this  paper  has  been  developed  to  improve 
the  speed  and  accuracy  with  which  the  minima  of  functions  ran  be  evaluated 
numerically.  In  addition,  a  matrix  characterizing  the  behavior  of  the  func¬ 
tion  in  the  neighborhood  of  the  minimum  is  determined  in  the  process. 
Linear  constraints  can  be  imposed  upon  the  variables  by  suitable  choice  of 
initial  conditions,  without  alteration  of  the  procedure. 


2.  NOTATION 

We  will  employ  the  summation  convention: 
N 

aM  bw  S  Y  by 
M  -  1 
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Ir.  dca*-ribing  the  iterative  procedure,  we  will  use  symbols  for  memory 
Ircatic-'s  rather  than  successive  values  of  .imimher;  e.g,,  we  would  write 
x  +■  3  — »•  x  instead  ol  x^  -r  3  =  x^4  j  In  this  notation,  the  sequence  of  oper 
aiiir.s  is  generally  relevant  The  f<  Ilowing  symbols  will  be  used. 


xH;  p  -  !  ,  X  the  set  of  N  independent  variables 

f  (x);  the  value  of  the  function  to  be  minimized  evaluated  at  the 
point  x- 

gu  (x):  the  di  riva fives  of  f  tx)  with  respect  to  x,u  evaluated. at  x: 

,  C’llX) 

gu  <51  = - ~ 

d  r-  x 


hu-;  d  non -r.egair. «.  symmetric  matrix  which  will  be  used  as  a 
metric  m  the  space  of  the  variables. 

A.  The  determinant  of  h1' ' 

£;  Z  times  fractional  accuracv  to  which  the  function  f  (&)  is 
to  be  minimized 

d:  a  limiting  value  for  who r  is  to  be  considered  as  a  "reason-, 
able"  minimum  value  of  the  function.  For  least  squares 
problems  d  can  be  set  equal  to  zero. 

an  integer  which  specifies  the  number  of  times  the  variables 
arc  to  be  changed  ir.  a  random  manner  to  test  the  reliability 
of  the  determination  of  the  minimum. 


K 


3.  GEOMETRICAL  INTERPRETATION 


It  is  convenient  to  use  geometrical  concepts  to  describe  the  mini¬ 
mization  procedure.  We  do  so  by  considering  the  variables  xM  to  be  the 
i ; :  rdir.ates  of  a  point  in  an  X-dimensional  linear  space.  As  shown  in 
Fig.  la,  the  set  of  x  for  which  f  (x)  is  constant  forms  an  N-l  dimensional 
s  .rface  in  this  space.  One  of  this  family  of  surfaces  passes  through  each  x, 
and  the  surface  about  a  point  is  characterized  by  the  gradient  of  the  function 
at  that  point: 


gM 


(*)  = 


These  N  components  of  the  gradient  can  in  turn  be  considered  as  the  coor¬ 
dinates  of  a  point  in  a  different  space,  as  shown  in  Fig.  lb.  As  long  as  f  (&) 
is  differentiable  at  all  points,  there  is  a  unique  point  g  in  the  gradient  space 
associated  with  each  point  x  in  the  position  space,  though  there  may  be 
more  than  one  x  with  the  same  g. 

r  [NOT  REPRODUCIBLE 


1II.A-3 


■  ■ 


INOT.  REPRODUCIBLE 


r 


> 


/ 


X' 


(a) 


(b) 


Fig.  1  .  Geometric.il  interpretation  of  xM  and  g^(x) 


In  the  neighborhood  of  any  one  point  A  the  second  derivatives  of 
f(x)  specify  a  linear  mapping  of  changes  in  position,  dx,  onto  changes  in 
gradient  dg,  in  accordance  with  the  equation 


dg.  = 


_JLl 

dxMdxV 


dx’ 


(3.1) 


The  vectors  dx  and  dg  will  be  in  the  same  direction  only  if  dx.is  an 
eigenvector  of  the  "Hessian  matrix: 

a2f  I 

CjX^C'X'1  ! 

If  the  ratios  among  the  corresponding  eigenvalues  are  large,  then  for  most 
dx  there  will  be  considerable  difference  in  the  directions  of  these  two 
vectors. 

All  iterative  gradient  methods,  of  which  this  is  one,  for  locating 
the  minima  of  functions  consist  of  calculating  j^for  various  x.  in  an  effort 
to  locate  those  values  of  x  for  which  g  =  0,  and  for  which  the  Hessjan 
matrix  is  positive  definite.  If  this  matrix  were  constant  and  explicitly 
known,  then  the  value  of  the  gradient  at  one  point  would  suffice  to  determine 
the  minimum.  In  that  case  the  change  desired  in  g  would  be  -j;,  so  we 
would  have 


dli 

c)x^ 


Ax17 


(3-2) 


from  wh’ch  we  could  obtain  Axv  by  multiplying  Eq.  (3.2)  by  the  inverse  of 
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the  matrix 


b\l 

cix^  dxv 


However,  in  most  situations  of  interest, 


II  cixM^x1*' 

is  not  constant,  nor  would  explicit  evaluation  at  points  that  might  be  far  from 
a  minimum  represent  the  best  expenditure  of  time. 


Instead,  an  initial  trial  value  is  assumed  for  the  matrix 


i)*f 


6x^dx1; 


This  matrix,  denoted  by  h^  ,  specifics  a  linear  mapping  of  all  changes  in 
the  gradient  onto  changes  ir  position.  It  is  Lo  be  symmetric  and  non-negative 
(positive  definite  if  there  are  no  constraints  on  the  variables?  After 
making  a  change  in  the  variable  x,  this  trial  value  is  improved  on  the  basis 

of  the  actual  relation  between  the  changes  in  g  and  x.  If  ^  ^ 


;s  con- 


11  <Jx^  ox 

stant, then,  after  N  iterations,  not  only  will  the  minimum  of  the  function  be 
determined,  but  also  the  final  value  of  hMV  will  equal 


o2f 


-l 


W( 


..  dxMdx*  u 

shall  subsequently  discuss  the  significance  of  this  matrix  in  specifying  the 
accuracy  to  which  the  variables  have  been  determined 


The  matrix  h  can  be  used  to  associate  a  squared  length  to  ony 
gradient,  defined  by  If  the  Hessian  matrix  were  constant  and  h^v 

were  its  inverse,  then  g^gv  would  be  the  amount  by  which  f(x)  would 

exceed  its  minimum  value  We  therefore  consider  as  specifying  a 
metric,  and  when  we  refer  to  the  lengths  of  vectors,  we  will  imply  their 
lengths  using  hr  as  the  metric  We  have  called  the  method  a  "variable 
metric"  method  to  reflect  the  fact  that  is  changed  after  each  iteration 


We  have  divided  the  procedure  into  five  parts  which  in  a  large  ex¬ 
tent  are  logically  distinct.  This  not  only  facilitates  the  presentation  and 
analysis  of  the  method,  but  it  is  convenient  in  programming  the  method  for 
machine  computation.  — 


*"■* INOT  'REPRODUCE 


4  READY.  CHART  1 


The  function  of  this  section  is  to  establish  a  direction  along  which 
to  search  for  a  relative  minimum,  and  to  box  off  an  interval  in  this  direc¬ 
tion  within  which  a  relative  minimum  is  located.  In  addition,  the  criterion 
for  terminating  the  iterative  procedure  is  evaluated. 

Those  operations  which  are  only  performed  at  the  beginning  of  the 
calculation  and  not  repeated  on  successive  iterations  have  been  included  in 
Chart  l  (page  7).  These  include  the  loading  of  input  data,  initial  print-outs, 
and  the  initial  calculation  of  the  function  and  its  gradient.  This  latter  cal¬ 
culation  is  treated  as  an  independent  subroutine,  which  may  on, its  initial  and 
final  calculations  include  some  operations  not  part  of  the  usual  iteration, 
such  as  loading  operations,  calculation  of  quantities  for  repeated  use,  special 
print-outs,  etc.  A  counter  recording  the  number  of  iterations  has  been  found 
to  be  a  convenience,  and  is  labeled I. 


i-J 
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INITIAL  CALL 


*Ootior>aHy  Printed  Slctenenfs 


The  iterative  part  of  the  compulation  uegins  with  'READY  l.a  The 
direction  of  the  first  step  is  chosen  by  using  the  metric  h  in  the  relation 

-h^gv— 

The  component  of  the  gradient  in  this  direction  is  evaluated  through  the 
relation 


g. 


(4.-2) 


From  Eqs.  (4.1.)  and  (4.2)  we  see  that  -gs  is  ‘be  squared  length  of  and 
hence  the  improvement  to  be  expected  in  the  function  is  -7gs.  The  positive 
definiteness  of  h^v  insures  that  g  -  is  negative,  so  that  the  step  is  in  a  direc¬ 
tion  which  (at  least  initially)  decreases  the  function.  If  its  decrease  is 
within  the  accuracy  deoired.  i  e  .  if  ga  '  €  >  0,  then  the  minimum  has  beet 
determined.  If  not, 'we  continue  with  the  procedure. 

In  a  first  effort  to  box  in  the  minimum,  we  take  a  step  which  is 

twice  the  size  that  would  locate  the  minimum  if  the  trial  hMvwere|  .  ^  ?jXv  ||  • 

However,  in  order  to  prevent  this  step  from  being  unreasonably  large  when 
the  trial  h^v  is  a  poor  estimate,  we  restrict  the  step  to  a  length  such  that 
(XsF)gn,  the  decrease  in  the  function  if  it  continued  to  decrease  linearly,  is 
not  greater  than  some  preassigred  maximum,  *(f-V).  We  then  change  x  by 


+  X 


(4.3) 


and  calculate  the  new  value  of  the  function  and  its  gradient  at  x  .  If  the 
projection  g+  =  g  +  of  the  new  gradient  in  the  direction  of  the  step  is 
positive,  or  if  the  new  value  of  the  function  P  is  greater  than  the  origina  ’ 
then  there  is  a  relative  minimum  along  the  direction  s  between  k  and  x  , 
and  we  proceed  to  "Aim"  w-here  wc  will  interpolate  its  position  However, 
if  neither  of  these  conditions  is  fulfilled,  the  function  has  decreased  and  is 
decreasing  at  the  point' x+,  and  we  infer  that  the  step  taken  was  too  small 
If  the  step  had  been  limited  by  the  preassigned  charge  in  the  action  d Re¬ 
double  d.  If  the  step  had  been  taken  on  the  hasp  of  1 i  ,  mo  1  y  ‘1  J° 
as  to  double  the  squared  length  of  s^,  leaving  the  length  of  a:i  pcrpendicu  ar 
vectors  unchanged  This  is  accomplished  by 


+  i  SM  sv 


(4.4) 


where  i  is  the  squared  length  of  sK  This  doubles  the  determinant  of  h^v. 
The  process  is  then  repeated,  starting  from  the  new  position. 


)  NOT  REPRODUCIBLE 


III.A-7 


5.  AIM:  CHART  L 


The  function  of  this  section  is  to  estimate  the  location  ot  the  rela¬ 
tive  minimum  within  the  interval  selected  by  "Ready.”  Also  a  ,  omr, arisen 
is  made  of  the  improvement  expected  b\  going  to  this  minimum  with  ihat 
from  a  step  perpendicular  to  this  direction. 
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CHART  Jt  AIM 


Inasmuch  as  the  interpolation  is  along  a  one-dimensional  interval, 
it  is  convenient  to  plot  the  function  along  this  direction  as  a  simple  graph 

(see  FiR.  i)- 

The  values  of  f  and  f+  of  the  function  at  points  x  and  x+  are  known; 
a,vd  so  are  its  slopes,  g*  and  g+.  at  these  two  points.  We  interpolate  for 
the  location  of  the  minimum  by  choosing  the  "smoothest”  curve  satisfying 
*e  boundary  conditions  at  x  and  x+.  namely,  the  curve  defined  as  the  one 

which  minimizes 
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Fig.  2 

Plot  of  f  (x)  along  a 
one -dimensional  interval. 


over  the  curve.  This  is  the  curve  formed  by  a  flat  spring  fitted  to  the 
known  ordinates  and  slopes  at  the  end  points,  provided  the  slope  is  small. 
The  resulting  curve  is  a  cubic,  and  its  slope  at  any  a  (0  )  is  given  by 


'  l  X  a  2 

Ss(  ;)  8g  ^  +  +  ^2”  (8s  +  8s  ^  ^z)> 


(5.1) 


where 


, .  lie'll + 


8s  *  8S 


T’ne  root  of  Eq.  (5.i)  that  corresponds  to  a  minimum  lies  between 
0  and  1  in  virtue  of  the  fact  that  gs  <  0  and  either  g£  >  0  or  z  <  gs  +  g+.  It 
can  be  expressed  as 


a  min  =  M1  -  a) 

where 

8s  f  Q  -  z 

8s  -  8S  +  2Q 

and 


Q  =  (**  -  8s8s)l/Z 


wot® 


i?.z) 


The  particular  form  of  Eq,  (5.2)  is  chosen  to  obtain  maximum  accuracy, 
which  might  otherwise  be  lost  in  taking  the  difference  of  nearly  equal 
quantities.  The  amount  by  which  the  minimum  in  f  is  expected  to  fall  be¬ 
low  f*  is  given  by 


P 

J(X-aX) 


5s (a)  =‘J  (gj  +  z  +  2Q)  a*  X 


(5.3) 
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The  anticipated  change  xs  new  compared  with  what  would  be  expected  from 
a  perpendicular  step.  On  the  basis  of  the  metric  hMv,  the  step  to  the  opti¬ 
mum  point  in  the  (N- 1  )-dimensional  surface  perpendicular  to  through 
x+^  is  given  by 


-h^g++-^-SM - -tM  (5.4) 

The  change  in  f  to  be  expected  from  this  step  is  jl‘J  g£.  Hence,  the 
decision  whether  to  interpolate  for  the  minimum  along  or  to  charge  x  by 
use  of  Eq.  (5.4)  is  made  by  comparing  =  t^  with  expression  (5.3) 

The  purpose  of  allowing  for  this  option  is  to  improve  the  speed  of 
convergence  when  the  function  is  not  quadratic.  Consider  the  situation  of 
Fig.  3.  The  optimum  point  between  x  and  x*  is  point  A.  However,  by  going 
to  point  B,  a  gr-cater  improvement  can  be  made  in  the  function.  When  the 
behavior  of  the  function  is  described  by  a  curving  valley,  this  option  is  oi 
particular  value,  for  it  ertables  successive  iterations  to  proceed  around 
the  curve  without  backtracking  to  the  local  minimum  along  each  step.  How¬ 
ever,  if  evaluation  of  the  function  at  this  new  position  does  not  give  a  hetter 
value  than  that  expected  from  the  interpolation,  then  the  interpolated  position 
is  used.  Should  the  new  position  be  better  as  expected,  it  is  then  desired  to 
modify  to  incorporate  the  new  information  obtained  about  the  function. 
The  full  step  taken  is  stored  at  s^,  and  its  squax-ed  length  is  the  sum  of  the 
squares  of  the  step  along  s  and  the  perpendicular  step,  i;e.,  sM=-gt+  +  X 2  £ , 
The  change  in  the  gradient  resulting  from  this  step  is  stored  at  g^s  and 
♦hese  quantities  are  used  in  the  section  "Dress"  in  a  manner  to  be  described. 


ax^+  (1  -  a)x+^—  (5.5) 

By  direct  use  of  the  xM  instead  of  the  aM  greater  accuracy  is  obtained  in 
the  event  that  a  is  small.  After  making  this  interpolation,  we  proceed  to 
••Fire." 
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6.  FIRE:  CHART  3 


The  purposes  of  this  section  are  to  evaluate  the  function  and  its 
gradient  at  the  interpolated  point  and  to  determine  if  the  lo<;al  minimum- 
has  been  sufficiently  well  located.  Tf  so,  then  the  rate  of  change  of  gra- 
dient  is  evaluated  (or,  more  accurately,  X  times  the  rate  of  change)  by 
interpolating  from  its  values  at  x,  x^,  and  at  the  interpolated  point. 

II  the  function  were  cubic,  then  f  at  the  interpolated  point  would 
be  a  minimum,  the  component  of  the  gradient  at  this  point  along  §_  would 
be  zero,  and  the  second  derivative  of  the  function  at  the  minimum  along 
the  line  would  be  2Q/X.  However,  as  the  function  will  generally  be  more 
complicated,  none  of  these  properties  of  f  and  its  derivatives  at  the  inter¬ 
polated  point  will  be  exactly  satisfied.  We  designate  the  actual  value  of  f 
and  its  gradient  at  the  interpolated  point  by  f  and  g.  .  The  component  of 
gp  along  s  is  g^  =  gs.  Should  f  be  greater  than  f  or  f+  by  a  significant 
amount  (>  c),  the  interpolation  is  not  considered  satisfactory  and  a  new  one 
is  made  within  that  part  of  the  original  interval  for  which  f  at  the  end 
point  is  smaller. 

From  the  values  of  the  gradient  gp,  g and  gp.  at  three  points  along 
a  line,  we  can  interpolate  to  obtain  its  rate  of  change  at  the  interpolated 
point.  With  a  quadratic  interpolation  for  the  gradient,  we 

(8M  -&)  +  tep  *  JM)  8 <61 ) 

where  ^  g^s  is  the  rate  of  change  of  the  gradient  at  the  interpolated  point 
The  component  of  gas  in  the  direction  of  s,  namely,  sf*  g^3  =  gss,  can  be 
expressed  as 


If  the  interpolated  point  were  a  minimum,  then  gs  ~  0  and  g8S  =  2Q. 

An  additional  criterion  imposed  upon  the  interpolation  is  that  the 
first  term  on  the  left  of  Eq.  (6.2)  be  smaller  in  magnitude  than  Q.  Among 
other  things,  this  insures  that  the  interpolated  value  for  the  second  deriva¬ 
tive  is  positive.  If  this  criterion  is  not  fulfilled,  no  interpolation  is  made, 
and  the  matrix  h^V  is  changed  in  a  less  sophisticated  manner 


l  NOT  REPRODUCIBLE 
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7.  DRESS.  CHART  4 


The  purpose  of  this  section  is  to  modify  the  metric  h^v  on  the  basis 
of  information  obtained  about  the  function  along  the  direction  jb.  The  new 
h;iV  is  to  have  the  property  that  (h^*')'  g^s  =  X  s^,  and  must  retain  the  infor¬ 
mation  which  the  preceding  iterations  had  given  about  the  function. 

If  the  vector  hH‘v  grs  =  t^  were  in  the  direction  of  sM,  then  it  would 
be  sufficient  to  add  to  hM*'  a  .matrix  proportional  to  sMsv.  If  tM  is  not  in  the 
direction  of  sM,  the  smallest  squared  length  for  the  dificrcnce  between  su 

X  1 

and  (hM v  t  asMsv)gvs  is  obtained  when  a  = -  -  -j  .  For  this  value  of  a, 

Ess  * 

the  squared  length  of  the  difference  is  t<j  -  -Esii  where  to  is  the  square 

& 

length  of  ji,  namely,  hM1'  d,j  dy.  When  this  quantity  is  sufficiently  small 
(<  e),  the  matrix  hMv  undergoes  the  change: 

ht*  t  (—  -  t  I  (7.1 ) 

The  corresponding  change  in  the  determinant  of  h^vis 


XK 

ESs 


A  A 


(7.2) 


When  the  vectors  tM  and  sM  are  not  sufficiently  colinear,  it  i9  necessary  to 
modify  hH^  by  a  matrix  of  rank  two  instead  of  one,  i.e., 


hMV.i^  + J_,M  sv. 
*0  8s  8 


,hMM 


(7.3) 


Then  the  change  in  the  determinant  of  is 


-Esa  A  — A  (7.4) 

lo 

After  the  matrix  is  changed,  the  iteration  is  complete;  after  printing  out 
whatever  infoimation  is  desired  about  this  part  of  the  calculation,  a  new 
iteration  is  begun.  This  is  repeated  until  the  function  is  minimized  to  with¬ 
in  the  accuracy  required 


8.  STUFF:  CHART  5 


The  purposes  of  this  section  are  to  test  how  well  the  function  has^ 
been  minimized  and  to  test  how  well  the  matrix  approximates  JJ  ^ 


at  the  minimum.  This  is  done  by  displacing  point  x  from  the  location  of  the 
minimum  in  a  random  direction. 


III.A-13 


16 


61  i 

63 

nmn 

— . — i 

u 

i  i 

RANDOM  NOS 

— V 

1_. 

l  “l 

f) 

— T—  — — 

I 

C«  L  1.  •„ 

^  •  **  *  A.*  •* 

AT  .** 

Rf  ADT  7 


FINAL 

PRINTOUT  I 

I  J 
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r“Al!  1  SV'F* 


ihe  dispfoccn-.cr.t  of  point  x  is  cho.cn  »bt  i  unit  length  in  terms 
of  fcPV  as  the  metric.  When  h-”  -  |!  >'  •  *“cl>  a  *“*  "iU  i""eas'-' 

f  by  half  the  square  of  the  length  of  the  step. 

If  the  directive  were  to  be  randomly  distributed,  then  iL  would  not 
be  satisfactuiy  to  choose  the  range  of  each  component  of  tM  independently; 
rather,  the  range  for  the  t„  should  be  such  that  h“v  t„  to  is  bounded  by 
preassigned  -.alues.  However,  this  refinement  has  not  been  incorporated 
into  the  charts  nor  the  computer  program.  The  length  of  the  step  has  been 
chosen  equal  to  one  so  that  the  function  should  increase  by  *  when  each 
random  step  is  taken. 

Significance  of  hJV: 

We  examine  a  least-squares  analysis  to  illustrate  how  the  initial 
trial  value  for  h^'  is  chosen,  and  what  its  final  value  signifies  In  this 
case,  the  function  to  be  minimized  will  be  chosen  to  beX  /*.  where  X  is 
the  statistical  measure  of  goodness  of  fit.  The  function  X  l  is  the  na  ura 
logarithm  of  the  relative  probability  for  having  obtained  the  observed  set 
of  data  as  a  function  of  the  variables  being  determined. 

The  matrix  h^  =  11  ^1^11  '  then  sPccifies  lhe  sPreads  and 
correlations  among  the  variables  by 

.  f  dNx  fxM.  <xM»(xv  -  <x’'  >)e_x2/2 

<  AxMAxv>  -  i - j~d Nx  e.xV2 


(8.1) 
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The  diagonal  elements  of  hf1*7  give  the  mean-square  uncertainty  for  each  of 
the  variables,  while  the  off-diagonal  elements  determine  the  correlations 
among  them.  The  full  significance  of  this  matrix  (the  error  matrix)  is  to 
be  found  in  various  works  on  statistics. (*)  It  enables  us  to  determine  the 
uncertainty  in  any  linear  function  of  the  variables,  for,  if  u  -  a^  x^  ,  then 


<u>  =  a^<  xM  > 

<  (Au  J2 >  -  (<x}1  xv>  -  <xh><xy  >) 

-  ,,  »i  '  '/**['  '.  (8-2a) 

If  u  is  a  more  general  function  of  then  if  in  a  Taylor  expansion  about  the 
value  of  x  derivatives  higher  than  first  can  be  ignored,  we  have 


<u  (x)  >  =  u  (<x>) 

<Aufa)>»  (  <&>  )5f5  «s»h" 


(6.  db) 


If  it  is  possible  to  estimate  the  accuracy  with  which  the  variables 
are  determined,  the  use  of  such  estimates  in  the  initial  trial  value  of  h 
will  speed  the  convergence  of  the  minimization  procedure.  Suppose,  for 
example,  that  to  fit  some  set  of  experimental  data,  it  is  estimated  that  the 
variables  xP  have  the  values: 


>:*  ^  3.6+  0  1 
x*  =  28.0  t  2 

x3  =  104±102  (83' 

Then,  the  initial  values  for  xM  and  hMv  would  be 

=  (3.0 

(0.01 

°o 

If  this  estimate  is  even  correct  to  within  a  couple  of  orders  of  magnitude, 
the  number  of  iterations  required  to  locate  the  minimum  may  be  substan¬ 
tially  less  than  that  for  some  more  arbitrary  choice,  such  as  the  unit 
matrix. 

If  it  is  desired  to  impose  linear  constraints  on  the  variables,  this 
can  be  readily  done  by  starting  with  a  matrix  h^v  which  is  no  lortger  posi¬ 
tive  definite,  but  which  has  zero  eigenvalues.  For  the  constraints 


28.0  1  04) 

0  0  \ 

4  0 

o  ioy  (8-4) 
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x^  -•  a 

V  =  2.  etc.. 

the  matrix  h^v  must  be  chosen  so  that 
hW  av  =  0 
h.m'  &*-  ■-  0 


(6.6) 


and  the  starting  value  for  x/(  .oust  satisfy  Eq .  (8  5).  For  example,  if  x3  is 
to  be  held  constant,  all  elements  of  hd'1'  in  the  third  row  and  third  column 
are  set  equal  to  zero  and  x3  is  set  equal  to  the  constant  value. 


When  constraints  arc  imposed,  instead  of  setting  A  equal  t'  the  de¬ 
terminant  of  h^v  (-0).  it  is  set  equal  lu  the  product  of  the  non-zero  eigen¬ 
value  of  hu'  .  Then,  except  for  round-off  errors,  not  only  will  the  <  onditions 
(8.6)  be  preserved  in  subsequent  iterations,  but  also  A  will  continue  toequal 
the  product  of  non- zero  eigenvalues 


Though  A  is  not  used  in  the  calculations,  its  value  may  be  of  interest 
in  estimating  how  well  the  variables  have  beoi  determined,  sin'x-c-h^f1  gives 


the  sum  of  the  eigenvalues  of  hu*\  while  A  gives  their  product.  T1  »•  square 
root  of  each  of  these  eigenvalues  is  equal  to  one  of  the  principal  st-tniaxes 
of  the  ellipse  formed  by  all  x_  for  which  f  (x)  exceeds  its  minimum  value  by 


9.  CONCLUSION 

The  minimization  method  described  has  been  coded  for  the  IBM-  <04 
using  Fortran.  Experience  is  now  being  gathered  on  the  operation  of  the 
method  with  diverse  types  of  functions.  Parts  of  the  procedure,  not  incor¬ 
porating  all  of  the  provisions  described  here,  have  been  in  use  for  some 
time  in  least- squares  calculations  for  such  computations  as  the  analysis 
of  Tf-P  scattering  experiments,^)  for  the  analysis  of  delayed  neutron  ex¬ 
periments,  (?)  and  similar  computations.  Though  full  mathematical  analysis 
of  its  stability  and  convergence  has  not  been  made,  general  considerations 
and  numerical  experience  with  it  indicate  that  minima  of  functions  car  be 
generally  more  quickly  located  than  in  alternate  procedures.  The  ability  of 
the  metric,  h^,  to  accumulate  information  about  the  function  and  to  oompen 
sate  for  ill-conditioned  is  the  primary  reason  for  this  advantage- 

The  author  wishes  to  thank  Dr.  G  Perlqw  and  Dr.  M.  Peahkin  for 
valued  discussions  and  suggestions,  and  Mr.  K.  Hillstrom  for  carrying  out 
the  computer  programming  and  operation. 
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APPENDIX* 

If  we  have  the  gradient  of  the  function  at  a  point  in  the  neighborhood 

||  »  then,  neglecting 

terms  of  higher  order,  the  location,  of  the  minimum  would  be  given  in 
matrix  notation  by 

£  =  x  -  G'1  V  (1) 

In  the  method  to  be  described,  a  trial  matrix  is  used  for_G_“l  and  a  step 
determined  by  Eq.  (1  }  is  laken.  From  the  change  in  the  gradient  resulting 
from  this  step,  the  trial  value  is  improved  and  this  procedure  is  repeated. 
The  changes  made  in  the  trial  value  for  _G"1  are  restricted  to  keep  t.hc  hunt¬ 
ing  procedure  "reasonable1’  regardless  of  the  nature  of  the  function.  Let 
H[  be  the  trial  value  for  GT1.  Then  the  step  taken  will  be  to  the  point 

x+  =  x  -  H  V  (2) 

The  gradient  at  x+,  7+,  is  then  evaluated.  Let  D  =  v+  -  V  be  the  change  in 
the  gradient  as  a  result  of  the  step  S  -  x+  -  x  =  -H  V.  We  form  the  new 
trial  matrix  by 

Hpp  =  HM.,+a(H^)^(H7+)v  (3) 

The  constant  a  is  determined  by  the  following  two  conditions: 

1  .  The  ratio  of  the  determinant  of  H+  to  that  of  JH  should  be  between 
R_1  and  R,  where  R  is  a  preassigned  constant  greater  than  I. 
This  is  to  prevent  undue  changes  in  the  trial  matrix  and,  in 
particular,  if  H  is  positive  definite,  H+  will  be  positive  definite 
also. 

2.  The  non- negative  quantity 

A  =  D  II +  D  +  S  (H  +  }_1  S  -  2  S  •  D  (4) 

is  to  be  minimized.  This  quantity  vanishes  when  S  =  H.+  D.  The  a  which 
satisfies  these  requirements,  together  with  the  corresponding  A,  as  functions 
of  N  -  V+  Hv+  and  M  =V+  HV  ,  are  as  follows:  (8) 


of  a  minimum  together  with  GT1,  where 


G=  II 


\ 

♦The  following  method  is  a  description  of  a  simplified  method  embody¬ 
ing  some  of  the  ideas  of  the  procedure  presented  in  this  report. 
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Range  of  M  j*_  A. 

M<-N/(R-1)  i/(M-N)  0  ,, 

-N/(R  -  1 )  <  M  <  N/(R  +  1 )  (l/RN)  -  (l/N)  (N  -  M  +  MR)l/RN 
N/(R  +  1)  <M  <NR/(R  +  1)  (N  -  2M)/N(M  -  N)  4M(N-M)/N 
NR/(R  +  1 )  <  M  <  NR/ (R  -  1 )  (R/N)  -  (l/N)  (M  ♦  NR  -  MR)  /RN 

NR/(R-1)<M  .  l/(M-N)  o  (5) 

The  dependence  of  Aon  M  is  bell- shaped,  symmetric  about  a  maximum  at 
M  =  N jZ,  for  which  a  =  0  and  A  =  N. 

After  forming  the  new  trial  matrix  H4,  the  next  step  is  taken  in 
accordance  with  Eq.  (2J  and  the  process  repeated,  provided  thatN  J  HV 
is  greater  than  some  preassigned  c.  When  the  gradient  is  constant,  it  can 

be  written  as 

V  =_G  (x  -  %  ).  <6) 

If  u  is  an  eigenvector  of  HG  with  eigenvalue  one,  then  it  will  be  an  eiger- 
vector  of  H+G  with  eigenvalue  one  as  well,  since 

H+Gu  =  HQ  u  +  a  H7  4  (!7+  HGu) 

=  u  +  a  H  V  4  [  v  HG  (1  -  HG)  u] 

=  u  •  (7) 

Furthermore,  when  &  -  0, 

H+  G  5  =  H/  D  =  S  (®) 

so  that  S  becomes  another  such  eigenvector.  After  no  more  than  N  steps 
(for  which  A  =  0).  H  will  equal  G"1  and  the  following  step  will  be  to  the  exact 

minimum . 

The  entire  procedure  is  covariant  under  an  arbitrary  linear  coordi¬ 
nate  transformation.  Under  these  transformations  of  x,V  transforms  as  a 
covariant  vector,  G  transforms  as  a  covariant  tensor  of  2nd  rank,  anri_H 
transforms  as  a  contra  variant  tensor  of  2nd  rank,  lhe  intrinsic  character¬ 
istics  of  a  particular  hunting  calculation  are  determined  by  the  eigenvalues 
of  the  mixed  tensor  HG,  and  the  components  of  the  initial  value  of  (x  -  £) 
aldng  the  direction  olThe  corresponding  eigenvectors.  Since  successive 
steps  will  bring  HG  closer  to  unity,  convergence  will  be  rapidly  accelerating 
even  when  G  itself  is  irregular.  Constraints  of  the  form  b  •  x  =  c  can  be 
improved  by  using  an  initial  H  which  annuls  b,  i.e., 

H  •  b  =  0 

and  choosing  the  initial  vector  x  such  that  it  satisfies  b  •  x  =  c.  Then  all 
steps  taken  will  be  perpendicular  to  b  and  this  inner  product  will  be  con¬ 
served.  For  example,  if  it  is  desired  to  hold  one  component  of  x  constant, 
all  the  elements  of  H  corresponding  to  that  component  arc  initially  set  equal 

to  zero. 
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